NASA Contractor Report 166528 



A Computer Program for the Calculation 
of Three-Dimensional Transonic 
Nacelle/Inlet Flowfields 


4IiAS&-Cfi-166S2d) k COfiPUTEB PfiOGEAM FOR THE 1184-13159 

CALCt}LA3:iON or THH££-EIH£H510!iAL IEANS029IC 
HACE1.LE/I8LE3: FROMFIELBS final Report 

ii^ct£eed-Georgia Co. « Marietta.} 102 p Onclas 

BC K06/M AOi CSCl 0 lA G3/02 43195 


Joseph Vadyak 
Essom H, Atto 

Contract NAS2-11285 
November 1983 






NASA Contractor Report 166528 


A Computer Program for the Calculation 
of Three-Dimensional Transonic 
Nacelle/lnlet Flowfields 


Joseph Vadyak 
Essam H. Atto 

Lockheed-Georgia Company 
Marietta, Georgia 


Prepared for Ames Research Center 
under Contract NAS2-11285 


rviASA 

National Aeronautics and 
Space Administration 


Ames Research Center 

Moffett Field. California 94035 




TABLE OF CONTENTS 

SECTION p^Q£ 

I. INTRODUCTION 2 

II. COMPUTATIONAL METHOD 4 

1. INTRODUCTION 4 

2. GRID TOPOLOGY AND GENERATION 4 

3. GOVERNING EQUATIONS FOR THE INVISCID FLOW . 8 

4. NUMERICAL ALGORITHM . . 11 

4.1 BASIC FINITE-DIFFERENCE SCHEME 11 

4.2 BOUNDARY CONDITIONS 14 

4.3 INITIALIZATION 15 

4.4 NUMERICAL STABILITY 16 

5. NUMERICAL RESULTS 16 

III. SUBROUTINE DESCRIPTIONS 33 

1. INTRODUCTION 33 

2. SUBROUTINE DESCRIPTIONS FOR THE NGRIDA 

GRID GENERATION PROGRAM 33 

3. SUBROUTINE DESCRIPTIONS FOR THE NACELLE 

FLOW SOLUTION PROGRAM 35 

IV. INPUT PARAMETERS AND OUTPUT INTERPRETATION ... 39 

1. INTRODUCTION 39 

2. NGRIDA GRID GENERATION PROGRAM INPUT 

PARAMETERS 39 

2.1 NAMELIST LISTl 39 

2.2 NAMELIST LIST2 . 47 

2.3 NAMELIST LIST3 48 

2.4 NAMELIST LIST4 50 

2.5 PARAMETER STATEMENT SPECIFICATION ... 51 

2.6 FILE USAGE 51 

3. NGRIDA PROGRAM OUTPUT INTERPRETATION .... 52 

4. NACELLE FLOW SOLUTION PROGRAM INPUT 

PARAMETERS 52 

4.1 NAMELIST LISTl 53 

4.2 NAMELIST LIST2 58 

4.3 NAMELIST LIST3 60 

4.4 NAMELIST LIST4 62 

4.5 PARAMETER STATEMENT SPECIFICATION ... 63 

4.6 FILE USAGE 63 

5. NACELLE PROGRAM OUTPUT INTERPRETATION. ... 63 

V. SAMPLE CASES AND SUGGESTIONS FOR USAGE 65 

1. INTRODUCTION 65 

2. NGRIDA GRID GENERATION PROGRAM 

SAMPLE CASES 65 


i 



2.1 SAMPLE CASE NO. 1 65 

2.2 SAMPLE CASE NO. 2 67 

2.3 SAMPLE CASE NO. 3 67 

2.4 SAMPLE CASE NO. 4 67 

3. NACELLE FLOW ANALYSIS PROGRAM SAMPLE CASES . 78 

3.1 SAMPLE CASE NO. 1 78 

3.2 SAMPLE CASE NO. 2 82 

3.3 SAMPLE CASE NO. 3 82 

3.4 SAMPLE CASE NO. 4 91 

4. SUGGESTIONS FOR USAGE 91 

4.1 PARAMETER STATEMENT USAGE 91 

4.2 NGRIDA PROGRAM USAGE 95 

4.3 NACELLE PROGRAM USAGE 95 

REFERENCES 97 


ii 



A COMPUTER PROGRAM FOR THE CALCULATION OF THREE-DIMENSIONAL TRANSONIC 

NACELLE / INLET FLOWFIELDS 


Joseph Vadyak and Essam H. Atta 
Advanced Flight Sciences Department 
Lockheed— Georgia Company 
Marietta, Georgia 30063 


SUMMARY 


A highly efficient computer analysis has been developed for 
predicting transonic nacelle/ inlet flowfields. This algorithm can compute 
the three-dimensional transonic flowfield about axisymmetric (or 
asymmetric) nacelle/ inlet configurations at zero or nonzero incidence. The 
flowfield is determined by solving the full-potential equation in 
conservative form on a body-fitted curvilinear computational mesh. The 
difference equations are solved using the AF2 approximate factorization 
scheme . 

This report presents a discussion of the computational methods used 
to both generate the body-fitted curvilinear mesh and to obtain the 
inviscid flow solution. Computed results and correlations with existing 
methods and experiment are presented. Also presented are discussions on 
the organization of the grid generation (NGRIDA) computer program and the 
flow solution (NACELLE) computer program, descriptions of the respective 
subroutines, definitions of the required input parameters for both 
algorithms, a brief discussion on interpretation of the output, and sample 
cases to illustrate application of the analysis. 

The isolated nacelle program developed herein has also been combined 
with the NASA-Ames TWING transonic wing flow analysis program to produce a 
wing/nacelle multicomponent flow analysis algorithm. The theory and usage 
of the wing/nacelle multicomponent program is discussed in a separate 
report. 



SECTION I 
INTRODUCTION 


Accurate Isolated nacelle/ Inlet flowfleld predictions are required 
for nacelle contour optimization for both axisymmetric and asymmetric 
configurations at zero and nonzero incidence. The prediction methods must 
be capable of analyzing three-dimensional transonic flowfields with the 
inclusion of embedded shock waves. The computational analyses which are 
developed must not only be accurate but also efficient enough to be used 
repetitively in parametric design studies. 


The objective of this investigation was to develop an accurate and 
highly-ef f icient method for calculating the three-dimensional transonic 
flowfield about axisymmetric (or asymmetric) nacelle/ inlet configurations 
at zero or nonzero incidence. The solution is obtained by solving the 
full potential equation in conservation law form on a three-dimensional 
body-fitted curvilinear mesh which is numerically generated. The use of 
the conservative form of the equation ensures that mass continuity is 
satisfied when capturing embedded shock waves. The difference equations 
are solved using the AF2 approximate factorization algorithm^. The AF2 
algorithm has been applied to the computation of two-dimensional transonic 
airfoil flows by Holst^ and to the computation of three-dimensional 
transonic wing flows by Holst and Thomas . Increases in convergence speed 
by factors of 4 to 7 have been realized using the AF2 scheme instead of 
using the standard transonic relaxation scheme, 
successive-line-over-relaxation. 

In the present investigation, the numerically— generated body— fitted 
grid is determined using a separate Fortran computer program called NGRIDA. 
The NGRIDA program can generate the computational mesh for axisynmetric 
nacelle/ inlet configurations. It can be modified to be applied in a 
meridional plane sense to obtain the computational mesh for asymmetric 
nacelle/ inlet configurations. 

The inviscid flow solution is obtained using a separate Fortran 
computer program called NACELLE. The NACELLE program uses the 
computational grid generated by NGRIDA or an alternate program. The 
NACELLE code can determine the flowfield solution for axisymmetric or 
asymmetric nacelle/ inlet configurations operating at zero or nonzero 
incidence in a subsonic or transonic free-stream flow. 

The development of the Isolated nacelle algorithm represents part of 
an overall effort to develop a combined wing/nacelle flow analysis program. 
The approach taken in developing the multicomponent algorithm is to use a 
component-adaptive grid embedding scheme in which the global computational 
grid is composed of a series of overlapped component grids, where each 
component grid is optimized for a particular geometry such as the wing or 
nacelle. The AF2 algorithm is used to determine the full-potential 
equation solution on each grid with trivariate interpolation being used to 
transfer property information between the component grids. To implement 
this approach for a wing/nacelle configuration, the NASA-Ames TWING 
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transonic wing flow analysis program^ has been combined with the NACELLE 
transonic nacelle flow analysis program presented herein. Preliminary 
results employing this scheme are reported by Atta and Vadyak^. a report 
detailing the theory and usage of the combined multicomponent algorithm is 
given in Reference 5. 
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1 . INTRODUCTION 

The three-dimensional body-fitted computational grid is determined 
using the NGRIDA computer program which employs a numerical mesh generation 
technique. The full-potential equation solution is obtained on the mesh 
generated by the NGRIDA program using the NACELLE flow analysis program. 

The efficient AF2 approximate factorization scheme is used by the NACELLE 
algorithm in solving the system of finite-difference equations. 

In this section, the computational methods used to generate the 
curvilinear mesh and to obtain the inviscid flow solution are presented. 
Computed results and comparison with existing axisymmetric flow analysis 
methods and three-dimensional flow experimental data are presented to 
illustrate application of the analysis. 


2, GRID TOPOLOGY AND GENERATION 

The flow solution is determined on a three-dimensional body-fitted 
curvilinear computational mesh. The computational mesh is obtained by 
using two-dimensional numerical grid generation techniques for a series of 
meridional planes which are splayed circumferentially around the body (a 
meridional plane is a plane containing the longitudinal axis of the inlet). 

Figure 1 illustrates the C— type nacelle/ inlet grid topology for a 
semi-infinite nacelle extended in the downstream direction. The base 
Cartesian coordinates are denoted by x, y, and z. Figure la shows the 
meridional plane grid topology, whereas Figure lb illustrates the grid as 
viewed along the longitudinal axis of the inlet which is the x axis. The 
computational curvilinear coordinates are denoted by ^ , n , and C . The C 
coordinate is in the wraparound direction, initiates at the external 
outflow surface, and terminates at the compressor face outflow surface. 

The n coordinate is in the circumferential direction, initiates at the ® » 

0 meridional plane, and terminates at the 6 » 2 tt meridional plane. The C 
coordinate is in the radial direction, initiates at the outer computational 
boundary (or centerline boundary), and terminates at the body surface. The 6 
- 0 and 6 - 27 f meridional planes are coincident for isolated nacelle 
configurations. If a pylon were present, then the pylon geometry would be 
contained between these two meridional stations. The (x-y) coordinate 
plane is the symmetry plane for angle of attack cases. 

The surface grid points are clustered in the region of the nacelle 
(leading edge point of the nacelle). The clustering is achieved by 
using a geometric stretching function which is expressed in terms of the 
arc length measured along the body from the hilite point in a given 
meridional plane. The stretching factor is used to compute the ratio of 
arc lengths between successive pairs of surface grid points as given by the 
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ct - AS /AS , 
s n n-1 


( 1 ) 


where AS^^ is the arc length between points n and n+1, and As_., is the arc 
length between points n“l and n^ as shown in Figure 2* Independent values 
of may be used for locating surface grid points on the external and 
internal nacelle contours* The outer computational boundary point 
distribution is determined by using either an angular distribution function 
expressed in terms of an angle measured about the nacelle hillte, or by 
again using an arc length distribution along the outer boundary* 


Once the surface and outer boundary mesh point distributions have 
been determined, the interior mesh point locations are computed using the 
NASA-Ames GRAPE® algorithm. The GRAPE algorithm determines the interior 
field point coordinates by iteratively solving two-coupled Poisson 
equations. For axisymmetric geometries, this two-dimensional grid 
generation procedure is applied for only one meridional plane, and the grid 
point locations on the remaining meridional planes are found by using 
simple reflection techniques. For asymmetric geometries, the 
two-dimensional grid generation algorithm can be applied for each 
meridional plane. 


Let X and y denote the Cartesian axial and radial coordinates, 
respectively, for a selected meridional plane, and let ^ and C denote the 
wraparound and radial curvilinear coordinates, respectively, for that 
plane. The interior grid point (x, y) distribution for the meridional 
plane is obtained by solving the following Poisson equations: 



- 

+ Yx 

CC 


x^ + Qx^) 

C2) 


- 267^^ 

+ Yy = 


+ Qy^) 

(3) 


In equations (2) and (3), the parameters a, 6, Y» and are defined as 

a “ X- V (4) 

(5) 


2^2 
X ^ + y 


® * *C ’'c 

Y - 


( 6 ) 


and 


M 


*5 


(7) 


The parameters P and Q are user— determined forcing functions which control 
grid point spacing and orthogonality near the computational boundaries. 
Equations (2) and (3) are solved numerically in the GRAPE algorithm using a 
successive-line-over-relaxation scheme. The boundary grid point locations 
are supplied to the GRAPE algorithm and are used as Dirichlet boundary 
conditions. Initial estimates of the interior field point locations are 
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required to start the relaxation procedure. 

The QRAPE two-dimensional grid generation algorithm is contained 
within the N6RIDA program and is used solely to obtain the Interior grid 
point locations. The boundary point locations are computed in NGRIDA and 
are supplied to the GRAPE algorithm to be used as boundary conditions. 

After the interior field point coordinates have been determined by GRAPE, 
point reordering, and grid translation, scaling, and reflection are 
performed in NGRIDA. The detailed theory and operation of the GRAPE 
algorithm are discussed in Reference 6. 

Figure 3 illustrates the rectangular computational space domain which 
corresponds to the physical space domain shown in Figure 1. The physical 
space boundaries transform into boundaries in computational space. This 
allows for accurate and straightforward boundary condition implementation. 
The mesh points are equally spaced in computational space thereby 
permitting the use of standard finite-difference formulae in the flow 
solution analysis. 

Figure 4 illustrates a typical meridional plane grid. The inlet 
centerline, which is the x axis, represents a singularity in the 
three-dimensional grid mapping. The computational grid boundary is 
displaced a small distance away from the x axis. An extrapolation and 
averaging technique, described later, is used to obtain flow properties on 
the centerline. Also, it has been found to be computationally advantageous 
to adjust the normal spacing near the external outflow surface to make a 
uniform radial point distribution along that surface. This produces a more 
favorable cell aspect ratio in this region, and thereby enhances stability 
and convergence in the flow solution algorithm. 


3. GOVERNING EQUATIONS FOR THE INVISCID FLOW 

The inviscid flow gas dynamic model is based on the assumption of 
steady potential flow which requires that the flow be both irrotational and 
isehtropic. The governing equations for steady three-dimensional potential 
flow are given by 


(^) + (^) + (^) - 0 

C n C 


[1 - 


Y-1 

Y+1 


(U <t>_ + V (^+W <p )] 


1 

Y-1 


( 8 ) 


(9) 


where C , n , and C denote the system of curvilinear coordinates, U, V, and W 
are the contravarianf velocity components in the Ci H» and c directions, 
respectively, p is the density, J is the Jacobian of transformation from 
the Cartesian coordinate system (x, y, z) to the general curvilinear 
coordinate system (C, H, ^), ^ is the velocity potential function, and Y is 
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Figure 4. Meridional plane grid. 
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the specific heat ratio. The density and contravariant velocity components 
are normalized by the stagnation density and critical sonic speed, 
respectively. 

Equation (8) is the full-potential equation in strong conservation 
law form. It expresses mass continuity for steady three-dimensional flows. 
Equation (9) expresses entropy conservation and is used to compute the 
density given the velocity potential field. 

The contravariant velocity components can be expressed in terms of 
the Cartesian velocity components as 


U = 

X y 


V = un + vn 
X y 


W = uC + 

X y 



(10) 

wn 

z 

(11) 

wC 

2 

(12) 


where u, v, and w denote the nondiraensional velocity components along the 
X, y, and z Cartesian coordinate axes, respectively. The contravariant 
velocity components can be also expressed in terms of the velocity 
potential (p as 


U = + Aj*^ 


V - A^+^ + Aj?^ + Ag*^ 


w = + A,<p + A^<p 

5 C 6 n 3 c 


The metric parameters (i 
transformation J are given by 


1 to 6) and the Jacobian of 


2 

2 


2 
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z 

6 XX 
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'z^: 

C n C 
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n ^ + 

C n c -> C n ; 




X y z 
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-C n c - C n C 

y X z X z y 


(13) 

(14) 

(15) 


(16) 


(17) 
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The following metric relations are also required in defining the mapping 
from physical space to computational space: 


X n c 


n = 

X 


- 

C = J(y -z 

X n 

Vs' 

^ = J(x^z - 

y c 0 


n = 
y 


- x^z^) 

C = J ( X^ 2 - - 

y n C 

(18) 


x.y^) 

n 

2 

J(x^y. 

- X ^ y ^ ) 

C = J ( X , y - 

2 n 


The metric 

parameters are 

obtained 

numerically using standard 



second-order or fourth-order accurate finite-difference formulae to compute 
derivatives of the form x_, x y^, etc. Then using the metric 

relations given by equations (18), the inverse quantities ? , H ^ etc., are 
determined. Substitution of these values into equations (I^) and (17) 
yields the and Jacobian J. 


4. NUMERICAL ALGORITHM 

4.1 BASIC FINITE-DIFFERENCE SCHEME 


The present numerical algorithm is based on the finite-difference 
formulation used by Holst and Thomas^ in computing transonic wing flows. 
In this algorithm, the full-potential equation residual is given by 


’-h.j.k * (f> 


i+4,j 


-4- n V 

n J 


i, j+4,k 


• 4 - 

6 

4 



i»j 


(19) 


where ^ denotes the residual operator, and i, j, and k. denote the 

grid poinf "^indices in the C (wraparound), n (circumferential), and 4 (radial) 
directions, respectively. The magnitude of the residual L(f). . ^ approaches 
zero as convergence is attained. The operators 5 r ( ) > and 5 ; ( ) 

are first-order accurate backward finite-difference operators (applied at 
midpoints) for the H, and directions, respectively. The terms P, "P, 
and P are upwind-biased density coefficients given by expressions of the 
form 
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^ ^i+h,J,k ^i+4i+r,j,k (20) 

where r denotes an upwind point along the K direction, and v is an 
artificial viscosity coefficient. In equation (20), the physical density p 
is computed using equation (9) with central differences used for 

determining the derivatives of 0. The artificial viscosity coefficient 
is given by 


V = 0 • M, . , < 1 

i.J ,k 


V = c(m: 


i» j ,k 


-i) 


M. > 


( 21 ) 


where M is the local Mach number, and C is a user-specified constant. The 
artificial viscosity coefficient C typically ranges from 1.0 to 2.0, with 
the larger values producing greater upwinding. Expressions similar to 
equation (20) hold for p and p which effect upwinding in the n and ^ 
d.;.cections , respectively. 


The finite-difference equations are solved using the AF2 approximate 
factorization scheme which has proved to be significantly more efficient 
than successive-line-overrelaxation schemes^. The AF2 algorithm is written 
in a three-step form as; 


First-Step: 


(ct 


It ^ 
T- 6 A.6 ) 
\ n J n 


’ij 


acjL(j) 


n 


( 22 ) 


Second-Step: 


\+i 


i,j ,k+l 


(A, + 6_ - - 6^A f^ , = 

k C 5 a C i r i,j,k Si,j 

Third-Step: 

(a + 6 )C^ . . = . . 

; i,d,k i,j»k 


(23) 


(24) 


12 



ORiG?f^Ai fan;: 

OF POOR QUAUTV 

In equations (22) to (24), a is a time-like factorization parameter chosen 
to maintain stability and attain fast convergence, 6 ^ is a factor that 
controls the amount of dissipation required in regions of supersonic flow, co 
is a relaxation factor, n is the iteration number, 14 ). ^ v is the mass 
residual [defined by equation (19)], f and g are inteniediate functions 
which are obtained during the solution process, and ^ is the 

potential function correction given by 


,n 

’i,j,k “ ^i,j,k 



(25) 


The terms Aj, and A^ are defined by 


A. 

1 


(pA^/J) 



,k 


(26) 


A. 

J 


(PA2/J) 


n 

i, 


(27) 


A 


k 


(PA 3 /J) 


n 

i.j ,k-ii 


(28) 


In Steps 1 and 2 the g and f functions are obtained by solving a 
tridiagonal system of equations while in Step 3 the correction C^. . , 
obtained by solving a bidiagonal system of equations. 


is 


The iterative relaxation procedure can be viewed as a time-marching 
integration algorithm in pseudotime. Using this analogy, the factorization 
parameter a appearing in equations (22), (23), and (24) can be regarded as 
the reciprocal of the marching pseudotime step.- To eliminate all 
components of the error frequency spectrum, it is generally preferable to 
employ a variable time step which sequentially varies with iteration 
number. The small values of a are particularly effective in reducing the 
low frequency errors, while the large values of a. are effective in reducing 
the high frequency errors. 


The factorization parameter a is computed using a repeating sequence 
which is based on iteration number^. The ot sequence that is the most 
frequently employed is given by 


k-1 


° “h (k = 1 . 2 . 


,M) 


(29) 


where is the lower limit of ct, is the upper limit, and is the 
value of a for the kth element of the sequence. For all cases presented 
herein, the number of elements in the sequence M was equal to 8. The 
optimum values of ot ^ and are generally determined by numerical 
experiment. For a (67xl3xl3)-point grid, factorization parameter limits of 
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■ 0.175 and » 6.0 were used. The relaxation parameter co is typically 
equated to 1.0. The time-'like dissipation factor 6^ is typically set to 
0.1. The use of 6^ ensures diagonal dominance in the ^ -difference equation 
tridiagonal matrix. The artificial viscosity coefficient C generally 
varies between 1.0 and 2.0. 


4.2 BOUNDARY CONDITIONS 

The density and potential function on the outer boundary (see Figure 
1) are held at their initial free-stream values during the course of the 
iterative solution procedure. As mentioned before, the grid centerline 
boundary is offset a small distance from the centerline (x axis) to avoid 
the mapping singularity. To determine property values on the centerline, 
an extrapolation and averaging procedure is employed. For each meridional 
pl^rie, the potential function values at interior mesh points are used to 
linearly extrapolate for the potential values on the centerline. A 
circumferential averaging of the extrapolated values is then performed for 
each wraparound station on the centerline. The potential values for points 
on each meridional plane centerline boundary are then updated by 
interpolation using property information from the centerline and interior 
field points. The updating of the properties on the grid centerline 
boundary is performed for every iteration. The property values on the 
outer and centerline boundaries are then used as Dirichlet boundary 
conditions for the ensuing iteration* Slight improvements in convergence 
speed can be realized by smoothing property values in the region of point 
(^) (see Figure 1), which delineates the transition point between the outer 
and centerline boundaries. 

At the body surface, the contravariant velocity W in the 
^—coordinate direction is specified. The velocity W is identically zero 
for a purely inviscid calculation as this satisfies the flow tangency 
condition. To account for boundary layer displacement effects, an 
effective can be calculated from the computed boundary layer growth and 
can then be applied as a surface transpiration velocity boundary condition. 
Once Wg is specified, the term is calculated from 

• (»S - ( 30 ) 


where the derivatives and are found using second— order differencing 
in the constant boundary surface. 


At both the compressor face and external outflow surfaces, the 
contravariant velocity in the ^“coordinate direction is specified. 
Given the local Cartesian velocity components u, v, and w, and the metric 


quantities 


Cy, and Cz, 


U 


Q can be determined using equation (10). At 


the compressor face, the flow is assumed to be uniform and in the axial (x) 
direction. The local flow velocity and density are fixed by specification 
of the required engine mass flow rate. At the external outflow surface, it 
is assumed that the free-stream velocity components are recaptured. This 


assumption is good far from the body surface, but is only approximate as 
the surface is approached. Locating the external outflow surface far 
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enough downstream minimizes the effect of this approximation. Once U has 
been determined, the term is calculated from ° 

* '“o ■ (31) 


where the derivatives c|)^ and 0 are found using second-order differencing 
in the constant boundary surface. 


At the mesh corner points ( C* C 


jflin ^ ^ ^ ^la^c ^ and ( ^ ^ inax ^ ^ ^ 

Cmax^* boundary conditions are enforced since both U and W are fixed. 

In these cases, 0 is computed using second-order differencing, and the 
derivative 0 is ^und from 




Aj - Aj Aj 


(32) 


The term (J)^ is then calculated using equation (31). 

In calculating the residual values for points on the body and outflow 
surfaces, the surface contravariant velocity along the direction normal to 
the surface is used directly. 

Boundary conditions are also required for the intermediate functions 
f and g. For isolated nacelle geometries, periodic boundary conditions are 
used in the n —coordinate direction. In the ^—coordinate direction, the 
condition fr “ is imposed at both the compressor face and external 
outflow surfaces. 


4.3 INITIALIZATION 


Significant improvements in convergence speed have been realized by 
initializing the potential function field using different procedures for 
the external flow and the internal flow. The potential field 
initialization for the external flow is performed using free-stream 
velocity components. For the internal flow initialization, the local 
velocity is assumed to be axial and is computed using a Mach number which 
is determined from the implicit relation 


M (1 


-(y^-1) 


= (~f) (1 ^ ^ 


-(y+1) 

2(y-1) 


(33) 


where M and denote the local and compressor face Mach numbers, 
respectively, and A and A^^ denote the local and compressor face flow 
areas, respectively. Equation (33) was obtained using one-dimensional gas 
dynamic formulae and simply expresses mass continuity for the captured 
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stream tube. To determine the internal potential field, a trapezoidal rule 
integration is used. The compressor face Mach number M . is fixed by 
specification of the engine mass flow rate. Generally, M ^ is directly 
entered in the program or is determined by specification of the inlet 
capture ratio. The free-stream Mach number M^j^^ and the compressor face Mach 
number M^^ can differ significantly at cruise conditions. 


4,4 NUMERICAL STABILITY 


During the course of program development, it was found necessary to 
impose a velocity potential function under-relaxation in the iterative 
process when calculating the solution at points near the centerline. The 
centerline represents singularity in the three-dimensional grid mapping. 
Without the under-relaxation, somewhat large potential function corrections 
would occur at mesh points near the centerline. 


A von Neumann linear stability analysis was conducted for the 
three-dimensional ,\F2 algorithm using the conservative form of the 
full-potential equation expressed in curvilinear coordinates. Simplified 
and linearized forms of the residual and AF2 correction operators were 
obtained by assuming that flow was subsonic (p ~ I), that there was a slow 
spatial variation of the metrics (A^^- 0. etc.), and that the grid was 
approximately orthogonal. Expressions were derived for the amplification 
factor in terms of Aj^ , A 2 , Aj, J, a, and w. The analysis indicated 
unconditional linear stability if 0 <"2 and ot > 0. When numerical values for 

the metric terms were substituted into the amplification factor 
expressions, low amplification factor values were obtained for the external 
flow mesh points and for the internal flow mesh points near the body 
surface. Higher amplification factor values were obtained near the 
centerline. Fastest convergence is attained if the amplification factor 
magnitude is near zero. 

To compensate for this effect, a simple potential function correction 
under-relaxation scheme was incorporated for mesh points close to the 
centerline. The potential function updating is thereby performed using 


n + 
i, j .k 


41. . , + ac. , , 

i» J 1, j ,k 


(34) 


where a ^1.0 at mesh points close to the centerline and ci = i.o elsewhere 
(standard algorithm). Generally, a is varied with the radial distance from 
the centerline. For mesh points close to the centerline, o is small (i.e., az 
.2), and as the radial distance is increased so is cr increased, finally 
reaching the value of 1.0 at a distance sufficiently far from the 
centerline. 


5. NUMERICAL RESULTS 

Selected numerical results are now presented to illustrate 
application of the analysis. Both axisymmetric and asymmetric 
nacelle/inlet configurations operating at zero and nonzero incidence are 


16 



considered. 


ORfGf^AL PASS [5^ 

OF POOR QUALHY 

The first set of computed results are for a recently-designed 
Lockheed-Georgia axisymmetric nacelle/ inlet configuration (GELACl inlet). 
Figure 5 illustrates the computed Mach number distributions (without 
viscous correction) for both the external and internal nacelle surfaces. 

The results shown are for a free-stream Mach number M^o of 0.8, an effective 
compressor face Mach number M^^ of 0.35, and an incidence ^ of 0 degrees. 
Plotted is the surface Mach number M versus the distance (AX) from the 
nacelle hilite. Since the geometry is axisymmetric and the angle of attack 
is zero, the flowfield will be axisymmetric. Also shown on Figure 5 are 
the results of the Jameson FLO-49 axisymmetric flow nacelle code^. The 
Jameson code is a two-dimensional full-potential finite-volume nacelle 
algorithm which uses successive-line-over-relaxation (SLOR) with multigrid 
convergence acceleration. It is limited to predicting the flow for 
axisymmetric nacelles at zero incidence. Referring to Figure 5, it can be 
seen that the two analyses agree very well. The Jameson program 
computation employed a grid with 128 wraparound stations and 32 radial 
stations. The three-dimensional nacelle program computation employed a 
grid with 67 wraparound stations, 13 circumferential stat ions ,• and 13 
radial stations. 

^-/hen an axisymmetric inlet is operated at incidence, cross flow 
develops and the flowfield becomes three-dimensional. Figure 6 illustrates 
the computed surface Mach number distribution when the GELACl inlet is 
operated at a = 3 (deg) while maintaining M = 0.8 and M ^ = 0.35. This 
figure illustrates the solution for the leeward ( 0= 0) , ^transverse [0 * 90 
(deg)]» and windward [0 » 180 (deg)] meridians. As the inlet is put at 
incidence, the stagnation point on the leeward meridian moves to increase 
the amount of flow expansion and subsequent peak Mach number and shock 
strength on the external surface, and to reduce the peak Mach number on the 
internal surface. The opposite effect occurs on the windward meridian. 

Figure 7 illustrates the computed surface Mach number distribution 
for a Lockheed-California axisymmetric inlet design (CALAC4 inlet) 
operating at ot * 0 (deg). Again, the free-stream and compressor face Mach 
numbers were maintained at 0.8 and 0.35, respectively. The results of the 
Jameson FLO-49 program are also presented on Figure 7 for the zero 
incidence case. The FLO-49 computation employed a (64 x I6)-point grid. 

Good agreement is obtained between the two analyses. 

Figures 8 through 12 show correlations between analysis and 
experimental data for the CALAC4 inlet. In Figure 8, the computed surface 
pressure coefficient C and wind tunnel data are compared for a case with M^ 
=» 0.6, a * 1.083 (deg), and an inlet capture ratio Q of 0.505. The capture*” 
ratio is defined as the ratio of free-steam capture area to the inlet 
hilite area. Figure 8 illustrates the results for both the external and 
internal surfaces for both the leeward and windward meridians. The 
agreement between the results of the numerical analysis and experiment is 
very good. The numerical results illustrated in Figure 3 were obtained 
using a computational grid with 67 wraparound stations, 25 circumferential 
stations, and 16 radial stations. The same case was executed using a 
(67X13X13)-point grid. The resulting surface solutions for the two grids 

ware almost identical except that the C peak on the windward meridian’s 

P 
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Figure 5. Solution for GELACl Inlet at M =0.8 and a • 0. 
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Figure 6. 


Solution for GELACl Inlet at = 0.8 and a =3 (deg). 
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Figure 7. Solution for CALAC4 Inlet at M =0.8 and a = 0. 
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Figure 8. Comparison of analysis and experiment for CALAC4 Inlet at M =0.6 
and a = 1.083 (deg) . “ 
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Figure 9. Comparison of analysis and experiment for CALAC4 Inlet at M ■ 0.601 
and a = 2.070 (deg) . ** 
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Figure 10. Comparison of analysis and experiment for CALAC4 Inlet at = 0.601 
and a = 3.066 (deg) . 
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CALAC4 INLET-LEEWARD MERIDIAN 



Fiffure 11 Comparison of analysis and experiment for CALAC4 Inlet at M = 0.851 
* * and a * 0.013 (deg). 
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Figure 12. Comparison of analysis and experiment for CALAC4 Inlet at M = 0.849 
and a = 1.074 (deg). “ 
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external surface vas better predicted using the finer grid. Figures 9 and 
10 compare analysis versus experiment maintaining the Mach number and 
capture ratio approximately constant but increasing the angle of attack, to a 
« 2.07 (deg) and 3.066 (deg), respectively. 


Figure 11 illustrates the computed surface pressure coefficient 
distribution and experimental data for the CALAC4 inlet operating at the 
conditions of - 0.851, a - 0.013 (deg), and Q - 0.2664. This represents 
a difficult test case because of the low capture ratio which requires that 
most of flow be spilled around the external cowling thereby generating a 
•large region of supersonic flow. The Mach number corresponding to the peak 
negative pressure coefficient is approximately 2*1. Nonetheless, excellent 
agreement is obtained with the data using a (67X13X13)-point grid. Figure 
12 illustrates the solution when the incidence has been increased to = 
1.074 (deg). These low capture ratio cases required increasing the 
^^tificial viscosity and time—like dissipation coefficients beyond their 
default values. This ensures numerical stability but decreases convergence 
speed. 


Figures 13, 14, and 15 illustrate the rapid convergence properties of 
the present approximate factorization algorithm. The variation of the 
number of supersonic points (NSP) with iteration number is plotted in 
Figure 13 for the flow field solution of the CALAC4 inlet operating at a 
free-stream Mach number of 0.8, a compressor face Mach number of 0.35, and ^ 
= 3 (deg) angle of attack. For this case, convergence, which corresponds 
to a three-order of magnitude reduction in maximum residual, was attained 
in 139 iterations. The value of NSP Increases very quickly with Iteration 
number. After 25 iterations, the NSP level has already achieved 90% of its 
final value. Figure 14 illustrates error reduction as a function of 
iteration number again for the 0^-3 (deg) incidence case of Figure 13. On 
the ordinate axis is plotted the log of the ratio of the absolute value of 
the maximum residual on the current iteration to that on the first 
iteration. (The results were plotted every 20 iterations). At 
^ppi^oxima tely 60 iterations, the maximum residual has been reduced by two 
orders of magnitude. Figure 15 presents the evolution of the surface 
solution with iteration number for the flow field computation of the GELACl 
inlet operating at - 0.8, M^^ - 0.35, and “ - 1 (deg). Shown in this 
figure are the computed solutions at 50, 80, and 158 iterations for the 
leeward and transverse meridians. The solution at 158 iterations 
represents the final converged solution at which point the maximum residual 
has been reduced by three orders of magnitude. After 50 iterations, a good 
approximation to the fully converged solution has been attained. The 
difference between the fully converged solution and that at 80 iterations 
is very minor. A characteristic of the AF2 algorithm is the rapid 
development of the surface Mach number and pressure distributions, as 
evidenced by Figure 15. 

Flow computations have also been performed for an asymmetric 
nacelle/ inlet configuration. The configuration under study is a recent 
Lockheed drooped-inlet design in which the inlet contour has 
circumferential variation both in section shape and in length from the 
hilite point to the compressor face. The front of the inlet is tilted 
downward with respect to the engine centerline for the purpose of aligning 
the inlet with the local flow direction underneath the wing. This reduces 
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Figure 15. Evolution of surface solution with iteration number. 
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cruise drag while maintaining the engine thrust vector in the proper 
orientation* Figure 16 illustrates the nacelle/inlet contours for the top 
(6 ■ 0) and bottom [6 - 180 (deg)] circumferential stations. The nacelle 
has one plane of geometric symmetry, and thereby the computational grid is 
generated for one half of the computational domain and then reflected about 
the symmetry axis to obtain the other half. 

For an asymmetric nacelle, the flow will be three-dimensional even at 
zero incidence. Computed surface Mach number contours for the top ( 6*0) 
and bottom [0 * 180 (deg)] meridians of the asymmetric nacelle are 

presented in Figure 17. These results are for a free-stream Mach number Moo 
of 0,8, a compressor face Mach number M^, of 0.35, and zero incidence ( a - 
0). Asymmetry exists in the computed solution with the lower external 
surface producing a higher local Mach number than the top meridian at this 
angle of attack. 

Typical computation times for three-dimensional flow solutions using 
the nacelle/inlet algorithm are 15-25 minutes for subcritical cases and 
35-50 minutes for supercritical cases on a VAX-11/780 super-mini-computer 
using a (67 x 13 x 13)-point computational grid. The algorithm has also 
been executed on the NASA-Ames CRAY-1 Class VI vector processor. Typical 
computation times for three-dimensional supercritical flow solutions using 
the CRAY-1 are on the order of 90 seconds of CPU time. No special effort 
has yet been made in vectorizing the program for use on the CRAY-1, 
Consequently, substantial improvements in the required execution time can 
be realized by restructuring the scalar code to benefit more fully from the 
CRAY vector-processing capabilities. It is estimated that by doing this 
the execution time could be reduced by a factor of 5 to 6 . 
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Figure 16. Asjnnmetric nacelle/inlet meridian contours. 
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Figure 17, Solution for asymmetric inlet at = 0.8 and a = 0. 
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1 . INTRODUCTION 

In this section, s brief description is given of the function of each 
subroutine in both the N6RIDA grid generation Fortran progrcun and the 
NACELLE flow solution Fortran program. This information supplements the 
information available in the form of comment statements within the 
respective programs. 


2. SUBROUTINE DESCRIPTIONS FOR THE NGRIDA GRID GENERATION PROGRAM 

NGRIDA . This program routine is the main control routine in the NGRIDA 
grid generation program. Routine NGRIDA first calls SUBROUTINE NINPUT for 
data input, parameter initialization, and printing of preliminary 
information. SUBROUTINE GENER is then called to control the actual 
numerical generation of the surface-fitted computational grid. 

♦ This subroutine is used to enter the input data, perform parameter 
initialization, and print preliminary information. The input data are 
entered by specification of the four NAMELISTS LISTl, LIST2, LIST3, and 
LIST4, All input parameters have specified default values, and a default 
axisymmetric nacelle/ inlet geometry data set has been loaded onto DATA 
statements within the program. At the end of SUBROUTINE NINPUT, the input 
namelist data are printed to initiate the computer output. SUBROUTINE 
NINPUT is called from routine NGRIDA. 

CUB ICS . This subroutine determines the coefficients for a set of cubic 
spline polynomials which are used to interpolate for the surface geometry 
given the axial positions and corresponding surface radii for N tabular 
data points. The use of cubic splines ensure continuity of surface radius, 
slope, and curvature for the Interpolated geometry. The spline 
coefficients are found in terms of the surface curvature values at each of 
the supplied tabular data points. A sparse system of simultaneous linear 
equations, which determines the curvature values, is solved by calling 
SUBROUTINE GELG. SUBROUTINE CUBICS is called from SUBROUTINE GENER. 

^NER. This subroutine contains the main control logic for numerically 
pnerating the curvilinear computational grid. The initial step in GENER 
is to generate the cubic spline polynomial coefficients for both the 
external and Internal nacelle surfaces by calling SUBROUTINE CUBICS. The 
cubic spline polynomials are used to interpolate for the surface geometry. 
At this stage, the arc length measured from the nacelle hilite is 
determined at each of the user-supplied tabular data points that define the 
geometry. From a set of user-supplied stretching factors which effect 
clustering around the hilite, the arc length distributions for the final 
surface grid point locations are determined. From the surface grid point 
arc length distributions, the axial and radial coordinates of the final 
surface grid points are found by interpolation using the cubic spline 
polynomials. At this stage, the outer computational boundary grid point 
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distribution is determined by one of two options. The first option 
determines the outer boundary point locations by using an angular 
distribution where the boundary points are located on rays spaced in equal 
angular increments around the nacelle hilite. The second option determines 
the outer boundary point locations by using an arc length distribution with 
the spacing between points being controlled by user-supplied stretching 
parameters. At this stage, a file is written for the GRAPE algorithm which 
determines the interior field point locations by iteratively solving two 
coupled Poisson equations. After the GRAPE algorithm has been executed, 
conversion of the data format and plotting of the finished grid is 
performed by calling SUBROUTINE CONVERT. At this stage, a reflection of 
the grid point coordinates is performed to determine the point 
distributions on the remaining meridional planes. Finally, certain 
transition point indices are determined and the grid point distribution is 
printed and also stored on a disk file for input to the flow solution 
program. SUBROUTINE GENER is called from routine NGRIDA. 

CONVERT . This subroutine is used to convert the output from the GRAPE 
algorithm into a format which is acceptable by the NACELLE flow solution 
algorithm. The GRAPE-generated grid point coordinates are entered through 
a binary read of TAPE 7. The corners on the leading edge of the grid are 
then rounded and point reordering is performed to make the grid generator 
output consistent with the input format required by the NACELLE program. 
Translation and scaling of the grid are then performed. If desired, 
plotting of the finished grid is performed. The plotting routines use the 
ISSCO-DISSPLA^ Fortran library of plot functions. If this subroutine 
library is not available, then the plotting section of the code will have 
to be modified to be compatible with the user’s particular installation. 
SUBROUTINE CONVERT is called from SUBROUTINE GENER. 

This subroutine is used for solving a system of simultaneous linear 
equations. The system is solved using Gaussian elimination with complete 
pivoting. This subroutine is IBM Library SUBROUTINE GELG. SUBROUTINE GELG 
is called from SUBROUTINE CUB ICS. 

The remaining subroutines in the NGRIDA grid generation program 
comprise the GRAPE algorithm. Cursory descriptions of the functions of 
these subroutines are presented herein. Detailed subroutine descriptions 
are given in Reference 6. 

This subroutine is the main control routine in the GRAPE algorithm. 
SUBROUTINE GRAPE is called from SUBROUTINE GENER after the boundary grid 
point locations have been determined. SUBROUTINE GRAPE determines the 
interior field point locations by calling in turn SUBROUTINES INPUT, 

INCHK, INNER, OUTER, SOLVE, and OUTPUT. The grid point locations are 
transmitted back to SUBROUTINE GENER through use of a scratch disk file. 

This GRAPE algorithm subroutine is used to check for smoothness in 
an input array by passing a parabola through the three nearest neighbors of 
a point and then observing the difference between the actual value and that 
interpolated from the parabola. 

CSPLIN . This GRAPE algorithm subroutine is used to perform cubic spline 
curve fitting. 
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This GRAPE algorithm subroutine Is used to initialize the grid 
coordinate and forcing function arrays prior to starting the Iterative 
relaxation procedure to determine the interior grid point locations. 


INCHK. This GRAPE algorithm subroutine is used to test for errors on 
selected input parameters for the GRAPE program. 


INNER. This GRAPE algorithm subroutine is used to locate grid points on 
the inner computational boundary. 


INPUT . This GRAPE algorithm subroutine is used to enter the input data 
that is required for execution of the GRAPE program. 

INTERP . This GRAPE algorithm subroutine is used to interpolate from the 
coarse grid solution to obtain initial conditions for the fine grid 
solution. 


OUTER . This GRAPE algorithm subroutine is used to locate grid points on 
the outer computational boundary. 

OUTPUT . This GRAPE algorithm subroutine is used to print the finished grid 
point coordinates. 

This GRAPE algorithm subroutine is used to solve the governing 
Poisson equations using a successive—line—over—relaxation iteration scheme. 

SOLVE . This GRAPE algorithm subroutine is used to control the coarse and 
fine grid solutions of the governing Poisson equations. 

This GRAPE algorithm subroutine is used to solve a tridiagonal 
system of simultaneous linear equations with fixed boundary conditions. 

This GRAPE algorithm subroutine is used to solve a tridiagonal 
system of simultaneous linear equations with periodic boundary conditions. 


3. SUBROUTINE DESCRIPTIONS FOR THE NACELLE FLOW SOLUTION PROGRAM 

NACELLE . This program routine is the main control routine in the NACELLE 
flow analysis program. Routine NACELLE first calls SUBROUTINE NINPUT for 
data input, initialization of certain parameters, and printing of 
preliminary information. SUBROUTINE NINIT is then called to perform 
selected parameter initialization, determination of outflow velocities, and 
potential field initialization. Finally, SUBROUTINE NSOLVE is called which 
applies the AF2 iteration scheme to obtain the Inviscid solution. 

NINPUT . This subroutine is used to enter the input data, perform 
initialization of selected parameters, and print preliminary information. 
The input data are entered by specification of the four NAMELISTS LISTl, 
LIST2 , LIST3, and LIST4. Most of the input parameters have specified 
default values. At the end of SUBROUTINE NINPUT, the input namelist data 
are printed to initiate the computer output. SUBROUTINE NINPUT is called 
from routine NACELLE. 


35 



CR-Q^?' AL PAQI: . • 

OF POOR QUALiTY 

NINIT * Thia subroutine is used to initialize selected parameters employed 
in the computation* SUBROUTINE NINIT first computes functions of the 
specific heat ratio and binomial expansion coefficients used in the 
computation of the physical density. SUBROUTINE METRIC is then called to 
determine the metric parameters used in the grid mapping from physical 
space to computational space. At this stage, the free-stream velocity 
magnitude and Cartesian velocity components are computed along with the 
appropriate outflow boundary condition parameters. Finally, the potential 
field Is initialized using one of two options. In the first option, the 
initialization is performed using velocity components based on the 
free-stream Mach number. In the second option, the external flow field 
points are initialized using free-stream velocity components, whereas the 
internal flow field points are initialized using velocities based on a 
local Mach number which itself is found from the area ratio and compressor 
face Mach number using one-dimensional gas dynamic formulae. SUBROUTINE 
NINIT is called from routine NACELLE. 

SECANT . This subroutine is used to determine the local internal flow Mach 
number using one-dimensional gas dynamic formulae which express mass 
continuity for the captured streamtube. A reference area and Mach number 
are supplied, and from the local flow area the local Mach number is 
determined. A secant iteration scheme is used to numerically determine the 
zero of the governing implicit relation. SUBROUTINE SECANT is called from 
SUBROUTINE NINIT and determines the compressor face Mach number from the 
user-supplied capture ratio. It also determines the local Mach number from 
the compressor face Mach number and the ratio between the local and 
compressor face areas. SUBROUTINE SECANT is called from SUBROUTINE NINIT. 

METRIC . This subroutine is used to compute the metric parameters for all 
points in the computational mesh. The metric coefficients define the 
three-dimensional grid mapping from physical space to computational space. 
The metric coefficients are computed using either second— order accurate or 
fourth-order accurate finite-difference formulae. The physical grid point 
locations are entered through a binary read of a disk file (TAPE 10). An 
option exists to calculate and print the finite-difference scheme 
amplification factor values as determined from a von Neumann linear 
stability analysis. SUBROUTINE METRIC is called from SUBROUTINE NINIT. 

NR0« This subroutine is used to calculate the physical density 
coefficients at the half points and the C-snd points in the computational 
mesh. From the current potential function field, the derivatives of the 
velocity potential in the 5, n» and C directions are computed using 
second-order accurate finite-difference formulae. The potential function 
derivatives are then used in a form of the isentropic energy equation to 
determine the density. At the body surface, the flow tangency condition is 
used in computing the density. At the internal and external outflow 
surfaces, the density is assumed known from the prescribed boundary 
conditions. At the external computational boundary, the density is assumed 
to be equal to the free-stream density. SUBROUTINE NRO is called from 
SUBROUTINES NSOLVE and NPRINT. 

• This subroutine is used to calculate the modified density 
coefficients at the C, 0, and C half-points in the computational mesh. The 
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modified density coefficients are computed by forming the appropriate 
averages of the physical density coefficients calculated by SUBROUTINE NRO. 
If the flow is locally supersonic, then upwindlng is effected by forming a 
weighted average of the physical density at the point in consideration and 
the density at a point in the upwind direction along the corresponding 
curvilinear coordinate* Upwindlng is always performed along the ^ 
wraparound coordinate. Upwindlng in the n circumferential direction and in 
the ^ radial direction are user-def eatable options. SUBROUTINE NROCO is 
called from SUBROUTINE NSOLVE. 

NRESID. This subroutine is used to apply the full-potential equation 
residual operator to the current velocity potential field to obtain the 
mass residuals at all points in the computational mesh except for those 
points on the outer computational boundary. The mass residual at a given 
point is computed by forming the sum of the finite-difference derivatives 
of the mass fluxes in the C.n , and ^ coordinate directions. Central 
difference expressions are used at interior field points while special 
procedures are used at the body surface and at the internal and external 
outflow surfaces. SUBROUTINE NRESID employs the modified density 
coefficients calculated by SUBROUTINE NROCO. SUBROUTINE NRESID is called 
from SUBROUTINE NSOLVE. 

NSOLVE . This subroutine is the main flow field integration control routine 
and is used to apply the AF2 approximate factorization scheme to determine 
the inviscid flow solution. The AF2 scheme is applied in an iterative 
manner with the iteration sequence being terminated by either achieving 
convergence or by reaching a maximum permissible number of iterations. Per 
iteration, a three-step process is used to obtain the potential function 
correction. The first two steps calculate intermediate potential function 
corrections arising from solving systems of factorized difference equations 
in the ^-wraparound and H -circumferential directions. The third step 
calculates the actual potential correction by solving a system of 
factorized difference equations for the C-radial coordinate. Steps 1 and 2 
require sweeping the computational mesh from the body surface outwards. 

The third step sweeps from the outer computational boundary towards the 
body surface. SUBROUTINE NSOLVE calls SUBROUTINES NRO, NROCO, and NRESID 
for calculation of the physical density coefficients, modified density 
coefficients, and mass residuals, respectively. SUBROUTINES NTRIC and 
NTRIP are called in the course of applying the AF2 scheme to solve systems 
of tridiagonal simultaneous linear equations with fixed and periodic 
boundary conditions, respectively. SUBROUTINE NPRINT is called to print 
the computed surface solution. If plots of the surface Mach number and/or 
pressure coefficient distributions are required, SUBROUTINE RESPLOT is then 
called by SUBROUTINE NSOLVE. SUBROUTINE NSOLVE is called by routine 
NACELLE. 


This subroutine is used to solve a system of simultaneous linear 
equations with a tridiagonal coefficient matrix and fixed boundary 
conditions. This routine is used in the AF2 iteration algorithm to obtain 
the potential function intermediate correction in the ^ -wraparound 
coordinate direction. SUBROUTINE NTRIC is called from SUBROUTINE NSOLVE. 

!■ This subroutine is used to solve a system of simultaneous linear 
equations with a tridiagonal coefficient matrix and periodic boundary 
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conditions. This routine is used in the AF2 iteration algorithm to obtain 
the potential function intermediate correction in the rj-icircumferential 
coordinate direction. SUBROUTINE NTRIP is called from SUBROUTINE NSOLVE. 

NPRINT . This subroutine is used to print the computed surface solution 
when convergence has been attained, or when the maximum permissible number 
of iterations has been reached, or throughout the computation at 
user-specified iteration counts. SUBROUTINE NRO is applied at the body 
surface to calculate the physical density coefficients at mesh mid-points 
on the surface. From the computed density distribution, the surface Mach 
number and pressure coefficient distributions are determined. SUBROUTINE 
NPRINT is called from SUBROUTINE NSOLVE. 

RES PLOT . This subroutine is used to plot the surface Mach number and 
pressure coefficient distributions at user-selected meridional stations. 
Options exist to defeat the plot option, or to plot just the surface Mach 
number distribution or the pressure coefficient distribution. The plotting 
routines use the ISSCO-DISSPLA^ Fortran library of plot functions. If this 
subroutine library is not available, then the plotting section of the code 
will have to be modified to be compatible with the user's particular 
installation. SUBROUTINE RESPLOT is called from SUBROUTINE NPRINT. 
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SECTION IV 

INPUT PARAMETERS AND OUTPUT INTERPRETATION 


1. INTRODUCTION 

In this section, the input parameters are defined for both the NGRIDA 
grid generation program and the NACELLE flow solution program* Where 
applicable, both the default value and typical value of each input 
parameter are given. Discussions on interpretation of the output, 

PARAMETER statement specification, and file usage for both programs are 
also provided. 

2. NGRIDA GRID GENERATION PROGRAM INPUT PARAMETERS 

The input data required for execution of the NGRIDA grid generation 
computer program are entered by namelist input. In all cases, the four 
NAMELISTS LISTl , LIST2, LIST3, and LIST4 are entered. All of the input 
parameters have assigned default values. A suitable grid can often be 
generated by retaining many of the input parameters at their default 
values . 


2.1 NAMELIST LISTl 


The input parameters specified in NAMELIST LISTl specify the number 
of mesh points, outer boundary shape, grid point distribution stretching 
factors, and print and plot options. 


The following three parameters specify the number of grid points. 

IMAX A positive integer denoting the number of ^-wraparound stations 

in the computational grid, as shown in Figure 18. The specified 
value of IMAX must be less than or equal to 140, and should be of 
the form (3nri-l) if grid sequencing is to be used in the GRAPE 
algorithm, where m is an integer. A default value of 67 is 
specified for IMAX. 

JMAX A positive integer denoting the number of n-circumferential 

stations in the computational grid, as shown in Figure 19. 
Stations J*1 and J»JMAX are coincident with the (x--y)-plane . The 
remaining meridional stations are spaced in equal angular 
increments around the body. The specified value of JMAX must be 
odd. A default value of 13 is specified for JMAX. 

KMAX A positive integer denoting the number of C-radial stations in 

the computational grid, as shown in Figure 18. The specified 
value of KMAX must be less than or equal to 60, and should be of 
the form (3n+l) if grid sequencing is to be used in the GRAPE 
algorithm, where n is an integer. A default value of 13 is 
specified for KMAX. 


The next three parameters specify the outer computational boundary. 
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XLEFT A real variable denoting the x~coordinate value of the left-hand 
side of the outer computational boundary, as shown in Figure 20. 
XLEFT is specified in the units of the input geometry, and 
typically is negative. A default value of -30.0 is specified for 
XLEFT. 

RADOUT A positive real variable denoting the radius of the outer 
computational boundary, as shown in Figure 20. RADOUT is 
specified in the units of the input geometry. A default value of 
50.0 is specified for RADOUT. 

DELTA A positive real variable denoting the offset distance between the 
computational boundary and the x-axis, as shown in Figure 20. 
DELTA is specified in the units of the input geometry and must 
be zero. A value of approximately 5% of the compressor face 
radius is recommended for DELTA. A default value of 0.5 is 
specified for DELTA. 

The next four parameters control the grid point distribution on the 

body surface. 

lOUTER A positive integer denoting the number of grid points desired on 
the outer portion of the nacelle surface ranging from the hilite 
point to the external outflow surface, as shown in Figure 20. A 
default value of 38 is specified for lOUTER. 

IINNER A positive integer denoting the number of grid points desired on 
the inner portion of the nacelle surface ranging from the hilite 
point to the compressor face outflow surface, as shown in Figure 
20. The specified values for lOUTER and IINNER must satisfy the 
relation 


IMAX =* IINNER + lOUTER - 1 


(35) 


A default value of 30 is specified for IINNER. 

ALPHAO A positive real variable denoting the stretching factor used in 
determining the grid point distribution on the nacelle external 
surface. Entering ALPHAO*!. 0 produces a point distribution with 
uniform arc length between the points [see equation (1)]. 

Entering ALPHAO > 1.0 clusters points near the nacelle hilite. A 
default value of 1.1 is specified for ALPHAO. 

ALPHAI A positive real variable denoting the stretching factor used in 
determining the grid point distribution on the nacelle internal 
surface. Entering ALPHAI»1.0 produces a point distribution with 
uniform arc length between the points [see equation (1)]. 

Entering ALPHAI > 1.0 clusters points near the nacelle hilite. A 
default value of 1.1 is specified for ALPHAI. 


The next five parameters determine the outer boundary grid point 
distribution. Generally, the outer boundary points near both the external 
and internal outflow surfaces are found by equating the axial coordinate of 
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Figure 20. Specification of outer boundary and point distribution on surface. 
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a given boundary point to the axial coordinate of the corresponding body 
point (see Figura 21). The remaining boundary points are found by either 
using an angular distribution or an arc-length distribution. 


lAA 


IDD 


KOUTER 


If KOUTER^ 
entered . 

ALPHBO 


ALPHBI 


A positive integer denoting the wraparound coordinate station up 
to which the points on the computational boundary near the 
centerline have their axial coordinates equated to the respective 
values of the corresponding points on the Internal nacelle 
surface, as Illustrated in Figure 21. lAA is specified in terms 
of the GRAPE algorithm point-ordering scheme in which the 
wraparound coordinate initiates at the compressor face and 
terminates at the external outflow surface (see Figure 21). A 
default value of 12 is specified for lAA. 

A positive integer denoting the wraparound coordinate station 
after which the points on the outer computational boundary have 
their axial coordinates equated to the respective values of the 
corresponding points on the nacelle external surface, as shown in 
Figure 21. IDD is specified in terms of the GRAPE algorithm 
point-ordering scheme in which the wraparound coordinate 
initiates at the compressor face and terminates at the external 
outflow surface (see Figure 21). Note that IDD must satisfy the 
relation lAA < IDD < IMAX. A default value of 57 is specified 
for IDD. 

An integer denoting whether the outer boundary grid points 
between points lAA and IDD are to be determined by an angular 
distribution or by an arc-length distribution. If KOUTER=0 , the 
grid points are determined by an angular distribution with the 
grid points being spaced at equal angular increments around the 
nacelle hilite. If K0UTER=1 , the grid points are determined by 
an arc-length distribution with respective stretching factors 
being specified by the user. A default value of 0 is specified 
for KOUTER. 

1 is specified, then the following two parameters must be 


A positive real variable denoting the stretching factor used in 
determining the outer boundary grid point distribution for the 
portion of the boundary between points C and IDD, as shown in 
Figure 21. Entering ALPHB0=1.0 produces a point distribution 
with uniform arc length between the points [see equation (1)]. 
Entering ALPHBO > 1.0 clusters points near point C. A default 
value of 1.0 is specified for ALPHBO. 

A positive real variable denoting the stretching factor used in 
determining the outer boundary grid point distribution for the 
portion of the outer boundary between point C and lAA, as shown 
in Figure 21. Entering ALPHBI»1.0 produces a point distribution 
with uniform arc length between the points [see equation (1)]. 
Entering ALPHBI > 1,0 clusters points near point C. A default 
value of 1.0 is specified for ALPHBI. 



ORIQJWaL PAGE SS 
OF POOR QUALITY 


liADIUS IG=IOO-l IG=iDO 



Figure 21. Definition of points lAA and IDD. 
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The next three parameters control the grid spacing adjacent to the 
body surface in the direction normal to the surface. Generally, the normal 
distance DS, as shown in Figure 22, is computed by the relation 

DS - (RCF - DELTA)/(KMAX - 1) (36) 


where RCF is the compressor face radius. For isolated nacelle grids, it is 
generally preferable to adjust DS near the external outflow surface to make 
the radial point distribution uniform as shown in Figure 22. This point 
distribution produces a more favorable cell aspect ratio at the outflow 
surface, and thereby enhances stability and convergence in the flow 
solution algorithm. 

DSFACT A positive real variable used as a multiplier for DS which is 
defined by equation (36). A default value of 1.0 is specified 
for DSFACT. 

KRAMP An integer denoting whether or not the DS distribution is to be 
altered to produce a uniform radial point distribution at the 
external outflow surface. If KRAMP =1 , the DS distribution is 
altered to produce the uniform radial point distribution. If 
KRAMP=0, the DS distribution is not altered, and the grid points 
will be packed towards the nacelle surface. If an isolated 
nacelle grid is being generated, then the KRAMP“1 option is 
recommended. A default value of 1 is specified for KRAMP. 

IRAMP A positive integer denoting the wraparound station after which 

the DS distribution is altered to produce a uniform radial point 
distribution at the external outflow surface, as shown in Figure 
22. IRAMP must be entered only if KRAMP- 1 . IRAMP is specified 
in terms of the GRAPE algorithm point ordering scheme in which 
the wraparound coordinate initiates at the compressor face and 
terminates at the external outflow surface. Note that IRAMP must 
satisfy the relation IRAMP < IMAX. A default value of 57 is 
specified for IRAMP. 

The next three parameters specify the scaling and smoothing function 
parameters. 


A positive real variable denoting the scaling multiplier used in 
obtaining the final grid. Entering SCALE-1.0 produces a grid in 
terms of the original input units. Setting SCALE =1 .0/RCF , where 
RCF is the compressor face radius, produces a grid with a 
compressor face radius of 1.0. This option is recommended as it 
allows for more easily determining the optimum acceleration 
parameters in the flow solution code. A default value of 
0.11904762-1.0/8.4 is specified for SCALE. 

An integer denoting the number of times smoothing polynomials are 
applied in determining the outer boundary axial coordinate 
distribution in the vicinity of point IDD, and the DS 
distribution in the vicinity of point IRAMP if the KRAMP-1 option 
is specified. If KTIME-0 is entered, no smoothing is performed. 
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Poor quality 


The default and recoomended value for KTIME is 3. 

NDEL A positive integer denoting the number of points to the left of 

point lOD at %diich the smoothing is initiated for the outer 
boundary point axial locations if KTIME ^ 0 is specified. NDEL 
must be entered only if KTIME > 0 is specified. A default value 
of 3 is specified for NDEL. 

The next two parameters specify the print options. 

KDUMP An integer denoting whether or not detailed nacelle surface and 
outer boundary coordinate data are to be printed. If KDUMP-0, 
the data are not printed. If KDUMP- 1 , the printout is performed. 
A default value of 1 is specified for KDUMP. 

KPRINT A one-dimensional integer array consisting of JMAX elements. 

Each element of KPRINT specifies whether or not a meridional 
plane grid point distribution is to be printed. Specifying 
KPRINT(J)=1 (J-1 , . . . , JMAX) causes the grid point coordinates to 
be printed for the Jth meridional plane. If KPRINT(J)-0 is 
entered, the printout is skipped for the Jth plane. A default 
value of KPRINT(l)-! is specified, while the remaining elements 
of KPRINT are specified as 0. 

The next two parameters specify the plot options. 

KPLOT An integer denoting whether or not the J=1 meridional plane grid 
is to be plotted. Specifying KPLOT-1 causes the grid to be 
plotted for the J-1 meridional plane. If KPLOT-0 is entered, the 
plotting is not performed. A default value of KPLOT-1 is 
specified. 

KDEV An integer denoting the plotting device. KDEV must be specified 

only if KPLOT-1 is entered. If KDEV-1 , the plot device is the 
VERSATEC electrostatic plotter. If KDEV— 2, the plot device is 
the HP ink plotter. A default value of 1 is specified for KDEV. 
(Note that the plot device specification can be altered by 
modifying SUBROUTINE CONVERT). 

The next parameter specifies the tape number (disk file) on which the 
grid coordinate data are written. 

ITAPE An integer denoting the tape (file) number on which the grid 

coordinate data are written. ITAPE can be entered as either 10 
or 14. For the present code version, ITAPE should be retained at 
its default value of 10. 

2.2 NAMELIST LIST2 

The input parameters entered in NAMELIST LIST2 specify the 
nacelle/ inlet geometry which is described by tabular input. A default 
axisymmetric nacelle/ inlet geometry has been loaded onto DATA statements 
within the program. The default geometry data are for the Lockheed-Georgia 
GELACl axisymmetric nacelle/ inlet configuration. 
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The following three parameters specify the internal surface contour 
for axlsymmetrlc nacelles. 


NI A positive Integer denoting the number of tabular data points 

used in defining the Internal surface contour, as shown in Figure 
23. A default value of 60 is specified for NI. 

XIN A one-dimensional real variable array consisting of NI elements. 

XIN(I) (I-1,...,NI) denotes the axial coordinate of the Ith point 
used in defining the internal surface contour, as shown in Figure 
23, Successive elements of XIN must be monotonically increasing. 
XIN(L) denotes the axial station of the nacelle hilite, and 
XIN(NI) denotes the compressor face axial station. Default 
values for XIN(I) (I-l,...,60) are specified. 

RIN A one-dimensional real variable array consisting of NI elements. 

RIN(I) (I«1,...,NI) denotes the radius of the Ith point used in 
defining the internal surface contour, as shown in Figure 23. 
RIN(l) denotes the nacelle hilite radius, and RIN(NI) denotes the 
compressor face radius. Default values for RIN(I) (I»l,...,60) 
are specified. * 


The next three parameters specify the external surface contour for 
axisymmetric nacelles. 

A positive Integer denoting the number of tabular data points 
used in defining the external surface contour, as shown in Figure 
23. A default value of 100 is specified for NO. 

XOUT A one-dimensional real variable array consisting of NO elements. 

XOUT(I) (I«1,...,N0) denotes the axial coordinate of the Ith 
point used in defining the external surface contour, as shown in 
Figure 23. The successive elements of XOUT must be monotonically 
increasing. XOUT(l) denotes the nacelle hilite axial station, 
and XOUT(NO) denotes the axial station of the external outflow 
surface. The relation XOUT(l) - XIN(l) must be satisfied. 

Default values for XOUT(I) (I-l 100) are specified. 

ROUT A one-dimensional real variable array consisting of NO elements. 

ROUT(I) (I»1,...,N0) denotes the radius of the Ith point used in 
defining the external surface contour, as shown in Figure 23. 
ROUT(l) denotes the nacelle hilite radius, and ROUT(NO) denotes 
the nacelle radius at the external outflow surface. The relation 
ROUT(l)»RIN(l) must be satisfied. Default values for ROUT(I) 

(I“l , . . . ,100) are specified. 

2.3 NAMELIST LIST3 


The input parameters entered in NAMELIST LIST3 specify the convergence 
criteria, iteration limits, and print options used in the GRAPE algorithm. 
These parameters are given only cursory descriptions herein, but are more 
fully discussed in Reference 6. Generally, the program is executed by 
retaining all of the parameters in NAMELIST LIST3 at their default values. 
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Figure 23, Surface geometry specification. 
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NORDA A one-dlttensional integer array consisting of two elements which 
denote the convergence criteria for the coarse-mesh and the 
fine-mesh GRAPE solutions. NORDA(l) and N0RDA(2) denote the 
numbers of order of magnitude by which the maximum correction is 
to be reduced for the coarse-mesh and fine-mesh solutions, 
respectively. Default values of N0RDA(l)-6 and NORDA(2)-3 are 
specified. 

MAXITA A one-dimensional integer array consisting of two elements which 
denote the iteration limits for the coarse-mesh and fine-mesh 
GRAPE solutions. MAXITA(l) and MAXITA(2) denote the iteration 
limits for the coarse-mesh and fine-mesh solutions, respectively. 
The MAXITA array is used in conjunction with the NORDA array. 
Default values of MAXITA(1 )=400 and MAXITA(2 )»200 are specified. 

JPRT An integer denoting the print option to be used in the GRAPE 

algorithm. If JPRT < 0, no printing is performed except for 
error warning messages. If JPRT»1 , a detailed printout will be 
performed. A default value of -1 is specified for JPRT. 

2.4 NAMELIST LIST4 

The input parameters entered in NAMELIST LIST4 specify the grid 
control functions and convergence parameters used in the GRAPE algorithm. 
Again, only cursory descriptions are provided herein, with detailed 
definitions being available in Reference 6. Generally, the program is 
executed by retaining all of the parameters in NAMELIST LIST4 at their 
default values. 


AAAI 

BBBI 


CCCI 

DDDI 


OMEGA 


OMEGP 

OMEGQ 


Positive real variables which control the enforcement of grid 
orthogonality and normal mesh spacing at the nacelle surface. 

Small values (i.e., 0.2) cause these effects to be propagated far 

into the field, whereas larger values (i.e., 0.6) cause these 
effects to diminish more rapidly. A default value of 0.6 is 
specified for both AAAI and BBBI. 

Positive real variables which control the enforcement of grid 
orthogonality and normal mesh spacing at the outer boundary. 

Small values (i.e., 0.2) cause these effects to be propagated far 

into the field, whereas larger values (i.e., 0.6) cause these 
effects to diminish more rapidly. A default value of 0.6 is 
specified for both CCCI and DDDI. 

A positive real variable used in the GRAPE successive— line— 
over-relaxation scheme to obtain the grid point distribution. 
OMEGA must be in the range of 0.0 < OMEGA <2.0. A default value 
of 1.3 is specified for OMEGA. 

Real variables used as relaxation parameters in obtaining the 
body surface forcing functions that are used in the Poisson 
equations. The effects of controlling grid spacing and 
orthogonality at the nacelle surface can be eliminated by 
entering OMEGP » OMEGQ =0.0. OMEGP and OMEGQ must be in the 
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range of 0.0 to 2.0. A default value of 0.1 is specified for 
both OMEGP and OMEGQ. 


OMEGR 

OMEGS 


Real variables used as relaxation parameters in obtaining the 
outer boundary forcing functions that are used in the Poisson 
equations. The effects of controlling grid spacing and 
orthogonality at the outer boundary can be eliminated by entering 
OMEGR ■ OMEGS • 0.0. OMEGR and OMEGS must be in the range of 0.0 
to 2.0. A default value of 0.0 is specified for both OMEGR and 
OMEGS. 


PLIM Real variables used as limitation factors in solving for the 

QLIM Poisson equation forcing functions. The range for each of these 

RLIM parameters is 0.0 to 100.0. A default value of 0.5 is specified 

SLIM for PLIM, QLIM, RLIM, and SLIM. 

DSOBI A positive real variable denoting the normal distance to be 

imposed between the radial stations KMAX and KMAX-1 . DSOBI has 
units of the input geometry. A default value of 0.2 is specified 
for DSOBI. 

2.5 PARAMETER STATEMENT SPECIFICATION 


The NGRIDA grid generation program uses variable array dimensions. 
The respective array sizes are fixed by specification of the following 
dimension parameters on the PARAMETER statement: 


Parameter 


Allowed Values 

NX 

> 

IMAX 

NY 


1 

NYG 

> 

JMAX 

NZ 


KMAX 

NXO 


NO 

NXI 


NI 

NXON 

> 

lOUTER 

NXIN 

> 

IINNER 

NXO 2 


N0*N0 

NXI2 


NI*NI 

s will be fixed 

at the 

time of compilatj 


2.6 FILE USAGE 


The following files are used by the NGRIDA grid generation code. 

File No. Usage 

TAPE 5 input file 

TAPE 6 printed output file 

TAPE 10 file on which the finished 

grid coordinate data are 
stored (output) 

TAPES 1,A,7 scratch files (should be 

deleted after computation) 



3. NGRIDA PROGRAM OUTPUT INTERPRETATION 
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The initial portion of the computer printout for the NGRIDA program 
consists of the NAMELIST input data. Then if the KDUMP-1 option is 
specified » detailed nacelle surface and outer boundary coordinate data are 
listed. After this, the GRAPE algorithm output is printed. The finished 
grid point coordinates are then printed for selected meridional planes as 
specified by the KPRINT array. Finally, certain index parameters are 
printed (whose definitions are given in the flow solution algorithm user’s 
manual) . 


The output parameters for the GRAPE algorithm are discussed in 
Reference 6. The remaining output parameters are defined below. 


I 

J 

K 

X 

Y 

Z 

XOUT 
ROUT 
XI N 
RIN 

ARC -LENGTH 
XBO 
RBO 
XBI 
RBI 
XSURF 
RSURF 
XBOUND 
RBOUND 


^-wraparound point index 
n-circumferential point index 
C-radial point index 

x-coordinate of point in finished grid, (scaled units) 

y-coordinate of point in finished grid, (scaled units) 

z-coordinate of point in finished grid, (scaled units) 

x-coordinate of input point on external surface, 

(original units) 

radius of input point on external surface, 

(original units) 

x-coordinate of input point on internal surface, 
(original units) 

radius of input point on internal surface, 

(original units) 

arc length measured from hilite point, 

(original units) 

x-coordinate of redistributed point on external 
surface, (original units) 

radius of redistributed point on external surface, 
(original units) 

x-coordinate of redistributed point on internal 
surface, (original units) 

radius of redistributed point on internal surface, 
(original units) 


surface point x-coordinate used as input for GRAPE 
code, (original units) 

surface point radius used as input for GRAPE code, 
(original units) 

outer boundary point x-coordinate used as input for 
GRAPE code, (original units) 

outer boundary point radius used as input for GRAPE 
code, (original units) 


4. NACELLE FLOW SOLUTION PROGRAM INPUT PARAMETERS 


The input data required for execution of the NACELLE flow analysis 
computer program are entered by namelist input. In all cases, the four 
NAMELISTS LISTl, LIST2, LIST3, and LIST4 are entered. Most of the inout 
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parameters have assigned default values. In many cases, the flow solution 
can be obtained by retaining most of the input parameters at their default 
values. The coordinate data for the computational mesh are entered through 
a binary read of TAPE 10. 

4.1 NAMELIST LISTl 

The Input parameters entered in NAMELIST LISTl specify the inlet 
orientation, free-stream and compressor face Mach numbers, convergence 
criterion, iteration limit, and Initialization, print, and plot options. 

The following six parameters specify the inlet orientation, free-stream 
Mach number, compressor face Mach number, and specific heat ratio. 

PITCH A real variable denoting the angle, in degrees, subtended by the 
free-stream velocity vector and the projection of that vector on 
the (x-z)-plane, as shown in Figure 24. A default value of 1.0 
degree is specified for PITCH. 

YAW A real variable denoting the angle, in degrees, subtended by the 

x-axis and the projection of the free-stream velocity vector on 
the (x-z)-plane, as shown in Figure 24. A default value of 0.0 
degrees is specified for YAW. 

XMFS A real variable denoting the free-stream Mach number. A default 

value of 0.8 is specified for XMFS. 

XMCF A real variable denoting the effective Mach number at the 

compressor face. XMCF must be entered only if the inlet capture 
ratio (CRATIO) is not entered. A default value of 0.3516 is 
specified for XMCF. 

CRATIO A real variable denoting the inlet capture ratio, which is 
defined by 

CRATIO » AINF/AHL (37) 

where AINF is the cross-sectional area of the capture streamtube 
far upstream of the inlet, and AHL is the effective hilite area 
(see Figure 25). If CRATI0i<0.0, then the compressor face Mach 
number XMCF will be internally computed from the supplied value 
of CRATIO. If CRATIO =* 0.0, then XMCF must be entered directly. 

A default value of 0.0 is specified for CRATIO. 

GAMMA A real variable denoting the specific heat ratio. A default 
value of 1.4 is specified for GAMMA. 

The next three parameters specify the metric calculation and potential 
function initialization options. 

KORDER An integer which denotes whether second-order or fourth-order 
accurate finite-difference formulae are used in obtaining the 
nistric parameters. If K0RDER=0, second-order accurate formulae 
are used. If KORDER = 1, fourth-order accurate formulae are 
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° IXJ ®°* (PITCH) cos (YAW) 
''oo = I^M" (PITCH) 


Wqo ~ Voo cos (PITCH) sin (YAW) 



Figure 24. Definition of PITCH and YAW angles. 
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HILITE POINT 



NACELLE 

Figure 25, Definition of AINF and AHL areas. 
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KINIT 


XMULM 

The 
limit . 

GRIT 

ITMAX 

The 

KPRCOR 

KPRINT 

SCALE 


used. The default and recommended value for KORDER Is 0. 

An Integer denoting the potential function initialization option. 
If KINIT ■ 0, the potential function initialization is performed 
using free-*stream velocity components for each point in the 
computational mesh. If KINIT ■ 1, the potential function 
initialization is performed using free-stream velocity components 
for the external flow mesh points, and using velocity components 
calculated from one*~dimensional gas dynamic formulae for the 
internal flow mesh points. The local internal flow Mach number 
is calculated using the compressor face Mach number and the ratio 
between the local flow area and the compressor face area. 

Setting KINIT - 1 generally significantly improves convergence 
speed. The default and recommended value for KINIT is 1. 

A positive real variable used as a Mach number multiplier for the 
internal flow points if the KINIT = 1 initialization option is 
specified. XMULM must be entered only if KINIT = 1, and 
generally is in the range of 0.6 to 0.8. The default value for 
XMULM is 0.7. 

next two parameters specify the convergence criterion and iteration 


A real variable denoting the convergence criterion used for 
terminating the calculation. Convergence is attained if the 
following relation is satisfied 


n 

I ^ax 



< GRIT 


(38) 


where is the maximum residual on the first iteration, and 

^mM maximum residual on the nth iteration. A default 

value of 0.001 is specified for GRIT. 


A positive Integer denoting the maximum number of iterations 
permissible. A default value of 200 is specified for ITMAX. 


next three parameters control the program output. 

An integer denoting whether or not the supplied grid point 
coordinates are to be printed. If KPRCOR-0, the mesh coordinates 
are not printed. If KPRCOR*!, the coordinate data are printed. 

A default value of 0 is specified for KPRCOR. 

A positive integer denoting the number of iterations between 
which the surface solution is printed. If KPRINT=1, the surface 
solution will be printed for each iteration. If KPRINT-50, the 
solution will be printed every 50 iterations, etc. If KPRINT is 
equated to a large number (i.e., 1000), then the solution will be 
printed only if convergence has been attained or if the Iteration 
limit has been reached. A default value of 1000 is specified for 
KPRINT. 


A positive real variable used as a multiplier for the surface 
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coordinates when the solution Is printed* This parameter allows 
the original coordinate units to be recovered if scaling was 
performed in generating the computational grid. A default value 
of 8.4 Is specified for SCALE. 

The next five parameters control the program plot options. 

KPLOT An integer denoting whether or not plotting of the surface 
solution is to be performed. If KPLOT-0, no plotting is 
performed. If KPL0T»1, plotting will be performed. A default 
value of 1 is specified for KPLOT. 

If KPL0T»1 is entered, the following four parameters must be specified. 

An integer denoting the plot device. If KDEV»1 , the plot device 
is the VERSATEC electrostatic plotter. If KDEV=2, the plot 
device is the HP ink plotter. A default value of 1 is specified 
for KDEV . (Note that the plot device specification can be 
altered by modifying SUBROUTINE RESPLOT.) 

A one-dimensional integer array consisting of JMAX elements, 
where JMAX is the number of circumferential stations in the 
computational grid (see NAMELIST LIST2). KMACH(J) ( J=1 , . . . , JMAX) 
specifies whether or not the surface Mach number distribution is 
to be plotted for the Jth meridional plane. If KMACH(J)»1, the 
Mach number distribution is plotted for the Jth meridional plane. 
If KMACH(J)“0, the plotting is not performed. All elements of 
KMACH have 0 specified as a default value except for KMACH(l), 
KMACH(4), and KMACH(7) which each have a default value of 1. 

A positive real variable denoting the maximum Mach number to be 
used in defining the ordinate axis for the Mach number plotting. 
XMTOP must be entered only if an element of KMACH is specified as 
nonzero. A default value of 1.5 is specified for XMTOP. 

A one-dimensional Integer array consisting of JMAX elements. 
KCP(J) (J*l , . . . , JMAX) specifies whether or not the surface 
pressure coefficient distribution is to be plotted for the Jth 
meridional plane. If KCP(J)»1, the pressure coefficient 
distribution is plotted for the Jth meridional plane. If 
KCP(J)-0, the plotting is not performed. All elements of KCP 
have 0 specified as a default value except for KCP(l), KCP(4), 
and KCP(7) which each have a default value of 1. 

The following parameter controls storage of the completed solution for 
use in post-processing or ensuing computations. 

ITDUMP An integer denoting whether or not the potential function field 
is to be loaded onto a disk file at the end of the computation. 

If ITDUMP-0, the file write is not performed. If ITDUMP-1 , the 
potential function values are loaded onto TAPE 12 (see SUBROUTINE 
NPRINT for the format of the write statement) . A default value 
of 0 is specified for ITDUMP. 


KDEV 


KMACH 


XMTOP 


KCP 
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4.2 NAMELIST LIST2 


The Input parameters entered in NAMELIST LIST2 specify the number of 
points in the computational mesh and certain index values. 

The following three parameters specify the number of mesh points. 

These parameters have the same definitions as noted in the grid generation 
algorithm user's manual. 

IMAX A positive integer denoting the number of ^-wraparound stations 
in the computational grid. A default value of 67 is specified 
for IMAX. 


JMAX 


KMAX 


A positive integer denoting the number of H-circumferential 
stations in the computational grid. A default value of 13 is 
specified for JMAX. 

A positive integer denoting the number of C-radial stations in 
the computational grid. A default value of 13 is specified for 
KMAX. 


The next parameter specifies whether certain index parameters are to be 
user specified or are to be read from the grid generator supplied disk 
file. 


ITREAD An integer denoting whether or not certain index parameters 

(defined later) are to be user specified or are to be read from 
the grid generator output disk file. If ITREAD-0, the index 
values must be user supplied. If ITREAD-1 , the index values are 
read from the disk file created by the grid generation algorithm. 
A default value of 1 is specified for ITREAD. 

If ITREAD"! is specified, then no other input parameters must be 
entered in this namelist. If ITREAD*0, the following four parameters must 
be entered (this allows for an alternate grid generation algorithm to be 
used) * 


A one-dimensional integer array consisting of KMAX elements. 
ITRANl(K) (K"I , . . . ,KMAX) specifies the ^—wraparound station that 
is closest to the radius of the nacelle hilite for the Kth radial 
station, as shown in Figure 26. The ITRANl array is used in 
determining the upwind direction for the C— wraparound coordinate. 
No default values are specified for the elements of ITRANl. 


A positive Integer denoting the ^-wraparound coordinate station 
at which the outer computational boundary becomes parallel to the 
inlet centerline, as shown in Figure 26. The ITRAN2 index is 
used to switch the boundary point treatment from that of a 
free-stream boundary to a centerline boundary. No default value 
is specified for ITRAN2. 


A positive integer denoting the C-wraparound coordinate station 
at which the potential function under-relaxation (described 
later) is to be initiated. Generally, ITRAN3 is between ITRAN2 
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ITRAN2 


Figure 26. Definition of ITRANl and ITRAN2 points. 
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and IMAX. No default value is specified for ITRAN3. 

IBEG A positive integer denoting the C-wraparound coordinate station 
at which the internal flow potential function initialization is 
started if the KINIT-1 option is specified. No default value is 
specified for IBEG. 

4.3 NAMELIST LIST3 

The input parameters entered in NAMELIST LIST3 specify the acceleration 
parameters, relaxation factors, artificial viscosity and damping 
coefficients, under-relaxation parameters, and smoothing functions used in 
the flow solution algorithm. 

The next four parameters specify the acceleration constants used in the 
AF2 iteration scheme. Two options have been incorporated into the program 
for calculating the acceleration parameter sequence. The first option 
employs the sequence 


K-1 

\ (K=1,2,3,...,M) (39) 

where is the lower limit for is the upper limit for ct, k is the 

sequence element number, and M is the number of elements in the sequence. 

The second option employs the sequence 


= ot 


H - 


K-1 
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(9-K) (K=2,3, . . .,7,8) 


(40) 


in which case M equals 8. 

KALPHA An integer denoting the ot -sequence to be used in the computation. 
If KALPHA-1, the sequence given by equation (39) is used. If 
KALPHA-2, the sequence given by equation (40) is used. The 
default and recommended value for KALPHA is 1. 


ALPHAL A positive real variable denoting the ot lower limit ot . 

A default value of 0.175 is specified for ALPHAL. ^ 

ALFHAH A positive real variable denoting the ct upper limit ot 
A default value of 6.0 is specified for ALPHAH. 

MD A positive integer denoting the number of elements M of the 

sequence in equation (39). MD must be entered as 8 if the 
KALPHA-2 option is specified. A default value of 8 is specified 
for MD. 


The following parameter specifies the standard relaxation factor used 
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OMEGA A positive real variable denoting the relaxation factor. 

*^or stability, OMEGA must be in the range 0.0 OMEGA < 2.0. 
.’he default and recommended value for OMEGA is 1.0. 


The next four parameters control the upwinding of the density and the 
amount of artificial viscosity and dissipation used in stabilizing the 
computations in regions of supersonic flow. Note that upwinding is always 
applied in the ^-wraparound coordinate direction. 


NJVIS An integer denoting whether or not upwinding in regions of 
supersonic flow is to be applied in the n^circumf erential 
coordinate direction. If NJVIS=0, the upwinding is not applied. 
If NJVIS=1, upwinding of the density is applied in the 
n-direction. Upwinding in the n-direction is needed only in 
cases of high cross flow. A default value of 0 is specified for 
NJVIS. 


NKVIS 


CFACT 


BXIE 


An integer denoting whether or not upwinding in regions of 
supersonic flow is to be applied in the ; -radial coordinate 
direction. If NKVIS=0, the upwinding is not applied. If 
NKVIS=1, upwinding of the density is applied in the C -direction. 

A default value of 0 is specified for NKVIS. 

A positive real variable denoting the artificial viscosity 
coefficient used in performing upwinding of the density in 
regions of supersonic flow. CFACT generally takes the values 1.0 
CFACT ^ 2.0, with higher values of CFACT producing greater 
upwinding. High values of CFACT can decrease shock wave 
resolution. A default value of 1.2 is specified for CFACT. 

A positive real variable denoting the timelike dissipation 
coefficient which is used in the AF2 algorithm to produce 
diagonal dominance in the^ —difference equations and thereby 
enhance stability in regions of supersonic flow. Increasing BXIE 
increases the amount of dissipation. A default value of 0.1 is 
specified for BXIE. 


The next three parameters control the smoothing operations performed in 
the computations. 


KSMTH 


KTIME 


A positive Integer denoting the ^ -radial station index K up to 
which smoothing is performed for the computed density in the 
^"Wraparound coordinate direction in the vicinity of the ITRAN2 
^—station. The density smoothing generally improves convergence 
speed. A default value of 4 is specified for KSMTH. 

An integer denoting the number of times a smoothing polynomial is 
applied in smoothing the density values in the C-direction in the 
vicinity of the ITRAN2 ^-station for the --stations ranging from 
1 to KSMTH. If KTIME=*0 is entered, no smoothing is performed. 

The density smoothing generally improves convergence speed. A 
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default value of 3 is specified for KTIME. 

NDEL A positive integer denoting the number of ^-stations to the left 
and to the right of point ITRAN2 at which the density is 
smoothed. Smoothing is performed from point (ITRAN2-NDEL) to 
point (ITRAN2+NDEL) for the ^“stations ranging from 2 to KSMTH, 
and from point (ITRAN2+1) to point (ITRAN2+2*NDEL) for the 1st 
^;-station. A default value of 3 is specified for NDEL. 

The next two parameters specify the potential function under-relaxation 

factors. 

CORFAT A one-dimensional real variable array consisting of KMAX 

elements. CORFAT(K) (K»l , . . . ,KMAX) specifies the potential 
function under-relaxation factor for the Kth C-station. This 
under-relaxation is necessary to maintain stability and enhance 
convergence for mesh points near the inlet centerline. 
Under-relaxation of the potential function is performed only for 
the ^-stations which are in the range ITRAN3 I ^ IMAX. The 
elements of CORFAT generally range from a low value at K * 1 to a 
value of 1.0 (no under-relaxation) as K approaches KMAX. The 
following default values have been assigned to the elements of 
CORFAT: CORFAT(K) = 0.2, 0.25, 0.30, 0.35, 0.40, 0.50, 0.60, 

0.70, 0.80, 0.90, 1.0, 1.0, and 1.0 for K = 1 to 13, 
respectively. 

ITCOR A positive integer denoting the iteration number at which the 

potential function under-relaxation is terminated. For iteration 
numbers greater than ITCOR, the standard algorithm with no 
under— relaxation is applied. Generally, fastest convergence is 
obtained by using under-relaxation for all iterations. The 
default and recommended value for ITCOR is 10,000. 

4.4 NAMELIST LIST4 

The input parameters entered in NAMELIST LIST4 specify if debug output 

is to be printed. 

KDUMP A one-dimensional integer array consisting of nine elements. 

Each element of KDUMP specifies whether or not a particular 
subroutine is to have debug output printed. Specifying 
KDUMP(I)«1 (I«l,...,9) activates the debug output option for the 
corresponding subroutine. Specifying KDUMP(I) = 0 causes the 
debug printout to be skipped. The elements of KDUMP activate the 
debug printout for the following subroutines/computations and 
have the following default values. 


KDUMP(I) 

Activates Debug Output 

Default 


for (SUBROUTINE) 

Value 

KDUMP(l) 

metric calculations (METRIC) 

0 

KDUMP(2) 

initialization (NINIT) 

0 

KDUMP(3) 

not presently used 

0 

KDUMP (4) 

AF2 solution scheme (NSOLVE) 

0 


62 



OF POOR QUALITY 

KDUMP(5) physical density calculation (NRO) 0 
KDUMP(6) residual calculation (NRESID) 0 
KDUMP(7) modified density calculation (NROCO) 0 
KDUMP(8) not presently used 0 
KDUMP(9) amplification factors (METRIC) 0 


ITSTRT 


An Integer denoting the iteration number at which debug output is 
to be initiated* A default value of 1 is specified for ITSTRT* 


4.5 PARAMETER STATEMENT SPECIFICATION 


The NACELLE flow solution program uses variable array dimensions. The 
respective array sizes are fixed by specification of the following 
dimension parameters on the PARAMETER statement: 


The array dimensions 
4.6 FILE USAGE 


Parameter 


NX ^ 

NY ^ 

NZ ^ 

NXP > 


will be fixed at 


Allowed Values 

IMAX 
UMAX 
raiAx 
IMAX + I 

time of compilation. 


The following files are used by the NACELLE flow solution code. 


File 

No. 

Usage 

TAPE 

5 

input file 

TAPE 

6 

printed output file 

TAPE 

10 

input file for grid coordinate 
data (from grid generation 
code) 

TAPE 

12 

potential function field 
(output) 


5. NACELLE PROGRAM OUTPUT INTERPRETATION 

The initial portion of the computer printout for the NACELLE flow 
analysis program consists of the NAMELIST input data* Then the transition 
P^i^t ii^dices ITRANl^ ITRAN2, ITRAN3, and IBEG are printed* After this, a 
convergence history summary is printed* Finally, the computed surface 
solution is output* The surface solution is determined at points midway 
between the original grid points and at the grid end points* 

The output parameters are defined below. 

I ^-wraparound point index 

J n “Circumferential point index 

^ x-coordinate of surface solution point 

Y y-coordinate of surface solution point 

Z z-coordinate of surface solution point 
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DX 


RHO 

MACH 

CP 

NSP 

RMAX 

RAVG 

IRMAX 

JRMAX 

KRMAX 


axial (x) distance between surface solution 
point and hilite point for a given meridional 
plane 

density coefficient defined as the local 

density normalized by stagnation density 

Mach number 

pressure coefficient 

number of supersonic points 

maximum residual 

average residual 

^-location of maximum residual 

n”locatlon of maximum residual 

C-location of maximum residual 


ommp± FAae m 
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SAMPLE CASES AND SUGGESTIONS FOR USAGE 


1 . INTRODUCTION 

In this section, a number of sample cases are presented to illustrate 
application of both the NGRIDA grid generation program and the NACELLE flow 
analysis program. For each sample case, a discussion of the problem Is 
given, the required input data are presented, and selected portions of the 
computer output are listed. The input parameter discussions follow the 
order in which the input parameters are presented in Section IV. At the 
end of this section, suggestions are provided for use of the respective 
computer programs. 

2. NGRIDA GRID GENERATION PROGRAM SAMPLE CASES 

Four sample cases are presented in this section to illustrate 
application of the NGRIDA grid generation algorithm to determining the 
computational grid about axisymmetric nacelle/ inlet configurations. 

2.1 SAMPLE CASE NO. 1 

This sample case is concerned with the grid generation for the 
Lockheed-Georgia GELACl axisymmetric nacelle/ inlet configuration. This 
sample case represents the case being considered when all input parameters 
are retained at their default values. 

The data deck for Sample Case No. 1 is illustrated in Figure 27. All 
parameters in NAMELIST LISTl are retained at their default values. 
Consequently, a 67x13x13 (IMAX x JMAX x KMAX) computational mesh will be 
generated. The outer computational boundary will have an outer radius 
RADOUT of 50.0, a left-side coordinate XLEFT of -30.0, and an offset 
distance DELTA of 0.5. The surface grid point distribution will be 
determined by the stretching factors ALPHAO and ALPHAI, which are both 
retained at their default values of 1*1, and by the distribution parameters 
IOUTER=38 and IINNER»30. The outer boundary grid point distribution is 
determined by retaining lAA and IDD at 12 and 57, respectively, and since 
K0UTER“0 , an angular distribution will be used in locating the points 
between points lAA and IDD. Since KRAMP*1 , the DS distribution will be 
adjusted to produce a uniform radial point distribution at the external 
outflow surface with IRAMP being retained at its default value of 57. All 
smoothing and plot options retain their default values. Since KDUMP-1 and 
KPRINT(1)»1, the nacelle and outer boundary grid points are printed along 
with the finished grid coordinates for the J*1 meridional plane. The 
meridional plane grid for the J-1 circumferential station will be plotted 
since KPL0T=1 . The finished grid will be scaled by the factor 
SCALE “0 . 11904762 which is the reciprocal of the input compressor face 
radius. Since ITAPE=10, the grid point coordinate data will be loaded onto 
disk file TAPE 10 for use by the NACELLE flow analysis program. 

All input parameters in NAMELIST LIST2 retain their default values. 
Consequently, the default nacelle/ inlet geometry (Lockheed GELACl inlet) 
will be used which has 60 internal surface tabular data points and 100 
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$LIST1 $END 
$LIST2 $END 
$LIST3 $END 
$LIST4 $END 


Figure 27. Data deck for NGRIDA Sample Case No. 1. 
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external surface tabular data points. The XIN, RIN, XOUT, and ROUT arrays 
are specified using the default geometry. 

All input parameters in NAMELIST LISTS and NAMELIST LIST4, which 
specify the GRAPE algorithm parameters, are retained at their default 
values. 

A plot of the completed grid is illustrated in Figure 28. Selected 
portions of the computer output for Sample Case No. 1 are presented in 
Figure 29. 

2.2 SAMPLE CASE NO. 2 

This sample case is concerned with the generation of a computational 
grid for the Lockheed-California CALAC4 axisymmetric nacelle. 

The data deck for this sample case is illustrated in Figure 30. All 
input parameters in NAMELIST LISTl are retained at their default values 
except for ALPHAO and ALPHAI which are now set as ALPHA0=l.I25 and 
ALPHAI=1.110. These larger values for the stretching factors enhance mesh 
resolution near the hilite. The nacelle/ inlet geometry is entered in 
NAMELIST LIST2. For this case, 64 and 83 tabular points are used in 
defining the internal and external nacelle contours, respectively; hence, 
NI=64 and NO=83. The surface coordinate data are then entered by the XIN, 
RIN, XOUT, and ROUT arrays. All input parameters in NAMELIST LIST3 and 
NAMELIST LIST4 are left at their default values. 

Selected portions of the computer output for Sample Case No. 2 are 
presented in Figure 31. 

2.3 SAMPLE CASE NO. 3 

This sample case is concerned with the grid generation for the 
Lockheed-California CALAC5 axisymmetric nacelle. 

The data deck for this sample case is illustrated in Figure 32. All 
parameters in NAMELIST LISTI retain their default values except for the 
mesh stretching factors which are again set as ALPHAO=1.125 and 
ALPHAI-1.110. The surface geometry is again entered in NAMELIST LIST2 with 
this case using 64 and 79 points in defining the internal and external 
contours, respectively. All input parameters in NAMELISTS LIST3 and LIST4 
retain their default assignments. 

Selected portions of the computer output for Sample Case No. 3 are 
presented in Figure 33. 

2.4 SAMPLE CASE NO. 4 

This sample case is concerned with generating a computational mesh for 
the CALAC4 inlet which has better mesh resolution than the grids generated 
in the previous sample cases. 

The data deck for this sample case is illustrated in Figure 34. Ail 
input parameters in NAMELIST LISTl retain their default values except for 
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Figure 28. Grid for NGRIDA Sample Case No. 1. 
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0\ 

vO 


listing 

01 

INPUT DATA 

SLTSTi 

TMAX 


67. 

JMAX 

s 

13. 

KMAX 

= 

13. 

XLEFT 

s 

-30.00000 , 

RADOIIT 

s 

50.00000 . 

delta 

= 

0,5000000 , 

TOUTFR 

s 

38. 

TINNER 

s 

30, 

ALPhAI 

a 

1.100000 . 

ALPHAO 

s 

1.100000 . 

kouter 

a 

0. 

ALPHPI 

a 

1.000000 , 

ALPHPQ 

a 

1,000000 . 

SCALE 

a 

0.1190476 . 

ttapf 

a 

10. 

KFLOT 

a 

1. 

KOEV 

= 

1 . 

kd«imp 

a 

t. 

TAA 

a 

12. 

TDD 

a 

57. 

KPPINT 

a 

i, 24»0, 

KTTNiE 

a 

3. 

NDEL 

a 

5, 

DSFATT 

a 

1. 000000 , 

KRAHP 

a 

1. 

tramp 

a 

57 

lEND 

SLT5T2 

NO 

a 

100, 

XOUT 

s 

O.OOOOOOOE+00, 4 

2.0000000F 

-02, 2.5000000F- 

4.5000002E 

-02, 5.0000001E-I 

7.0000000E 

-02, 7,999999BF-i 


(^1 


>999F-03, 9,999999bE-03, 1 . 5n0Oo00K-02 , 

2,9999999F!-02, 3 . 5000000F-n2 , 3 . 9999999E-02 , 

5,b0O00O0R-02, 5 . 9999999F!-02 , 6 . 4999998 E-02 , 

9.0000004E-02, 0.1000000 , 0.1100000 . 


Figure 29. Output for NGRIDA Sample Case No. 1. 
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o 


Ix 

1 

J3 

1 

Ks 

1 

Xs 

O.R33333F+01 

Ys 0.b9b?38b:+01 

Zs 

O.OOOOOOE+00 

Ic 

2 

J« 

1 

Ks 

1 

Xs 

0.792466E'+01 

Ys 0.b9b238c:-»'01 

Zs 

0.000000C<f00 

IB 

3 

J» 

1 

Ks 

1 

Xs 

0.752527F:+Oi 

Ys 0,b9b238h>01 

Zs 

o.noooooE+oo 

Is 

4 

JB 

1 

Ks 

1 

Xs 

O.715054F+01 

Ys G.b9b2i8t+01 

Zs 

O.OOOOOOC^OO 

1b 

5 

J* 

1 

Ks 

1 

Xs 

0.660709F:+01 

Ys C.b9b238t: + 01 

Z = 

0.O00000E400 

Ib 

6 

J* 

1 

Ks 

1 

Xs 

0.fi49343E+01 

Ys 0.b9b238E-^0i 

Zs 

0,OOOOOOE<fOO 

Ix 

7 

d* 

1 

Ks 

1 

Xs 

0.620492F+01 

Ys 0.b9b238t-»-01 

Zs 

r.ooooooE^oo 

Ib 

R 

JB 

1 

Ks 

1 

Xs 

0.593510F+01 

Ys O.59b230E+Ol 

Zb 

O.OOOOOOE4>00 

Ib 

9 


1 

Ks 

1 

Xs 

0.567372P+01 

Ys 0.59523RE+01 

Zs 

O.OOOOOOE^OO 

Ix 

10 

d» 

1 

Kx 

1 

Xs 

0.540402F+01 

Ys 0.59b238E+01 

Zs 

0.OO0000E400 

Is 

11 

J3 

1 

Ks 

1 

Xs 

0.510106E+01 

Ys 0,b95?3RE+01 

Zs 

O.OOOOOOE400 

Ib 

12 

JB 

1 

Ks 

1 

Xs 

0.473751E+01 

Ys 0.595238E+01 

Zb 

O.O00000E400 

Is 

13 

d* 

1 

Ks 

1 

Xs 

0.430742E+01 

Ys 0.b95238E+01 

Zs 

O.OOOOOOE^OO 

Ic 

14 

JB 

1 

Ks 

1 

Xs 

0. 384002F-f01 

Ys 0.59b23RG+01 

Za 

0,ooOOOOE4>00 

Is 

15 

J3 

1 

Ks 

1 

Xs 

0.336708F+01 

Ys 0.b9523BE<f01 

Z® 

O.OOOOOOEtOO 

Is 

IR 

Js 

1 

Ks 

1 

Xs 

0.?«9572F4-01 

Ys 0.b9b23Pt+01 

Zs 

O.OOOOOOEtOO 

Is 

17 

U® 

1 

Ks 

1 

Xs 

O.242287Et0i 

Ys 0.5Ob238E+0l 

Zs 

O.OOOOOOE<fOO 

Is 

18 

JS 

1 

Ks 

1 

Xs 

0.192223E+01 

Ys 0.59523Rt:*Ol 

Zb 

O.OOOOOOEtOO 

Is 

19 

J= 

1 

Ks 

1 

Xs 

0.1 39033F401 

Ys 0.b9b23Re+01 

Zs 

O.OOOOOOE^OO 

Is 

20 

JS 

1 

Ks 

1 

Xs 

O.R56567E-^00 

Ys G.b90459E+0l 

Zs 

O.OOOOOOE<fOO 

Is 

21 

J® 

1 

Ks 

1 

Xs 

O.439945F+00 

Ys C.568599E+01 

Zs 

O.OOOOOOE4>00 

Is 

22 

JS 

i 

Ks 

1 

Xs 

0.83881 7F-01 

Ys 0.b?1326E+0i 

Zs 

O.OOOOOOE^OO 

Is 

23 

JS 

1 

Ks 

t 

Xs 

O.Oo^OOOEtOO 

Ys 0.454335E+U1 

Zs 

O.OOOOOOE+00 

Is 

24 

JS 

1 

Ks 

1 

Xs 

O.OuOOOOF+00 

Ys 0,392977E+0t 

Zs 

O.OOOOOOE+OO 

Ix 

25 

JS 

1 

Ks 

1 

Xs 

O.nooOOOE-fOO 

Ys 0.34050ftt4.01 

Zs 

0,Oo0000E4>00 

Is 

2fi 

Js 

1 

Ks 

1 

Xs 

0,OOOOOOF4>Oo 

Ys 0.294392E+01 

Z® 

O.OQOOOOE4>00 

Ts 

27 

JS 

1 

Ks 

1 

Xs 

O.POOOOOE+OO 

Ys G,252R83E+01 

Zs 

O.OOOOOOE+00 

Is 

28 

JS 

1 

Ks 

1 

Xs 

O.OOOOOOE+00 

Ys G.214705E+01 

Zs 

O.OOOOOOE<fOO 

Is 

29 

J* 

1 

Ks 

1 

Xs 

0,O00000F+O0 

Ys 0,17HR8flE+0t 

Zs 

O.OOOOOOE+OO 

Is 

30 

Jx 

1 

Ks 

1 

Xs 

0,00O000F+00 

Ys 0.144R57E401 

Zs 

O.OOOOOOEtOO 


O O 


*0 
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V3 r- 

O 

c 

> 
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m 






Figure 29. Continued. 
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$LIST1 ALPHAO-1.125, ALPHAI-1.11 $END 
$LIST2 NI-64, NO-83, 


XIN( 

D- 

0.0000 

,RIN( 

D- 

7.6822, 

XIN( 

2)- 

0.0016 

,RIN( 

2)- 

7.6498, 

XIN( 

3)- 

0.0065 

,RIN( 

3)- 

7.6173, 

XIN( 

4)* 

0.0147 

,RIN( 

4)- 

7.5849, 

XIN( 

5)- 

0.0261 

,RIN( 

5)- 

7.5525, 

XIN( 

6)- 

0.0410 

,RIN( 

6)- 

7.5200, 

XIN( 

7) = 

0.0593 

,RIN( 

7) = 

7.4876, 

XIN( 

8) = 

0.0811 

,RIN( 

8)- 

7.4551, 

XIN( 

9)- 

0.1066 

,RIN( 

9)- 

7.4227, 

XIN( 

10) = 

0.1359 

,RIN( 

10)- 

7.3903, 

XIN( 

11)« 

0.1693 

,RIN( 

11)- 

7.3578, 

XIN( 

12)- 

0.2068 

,RIN( 

12)- 

7.3254, 

XIN( 

13)- 

0.2488 

,RIN( 

13)- 

7.2929, 

XIN( 

14)- 

0.2957 

,RIN( 

14)- 

7.2605, 

XIN( 

15)- 

0.3477 

,RIN( 

15)- 

7.2281, 

XIN( 

16)- 

0.4055 

,RIN( 

16)- 

7.1956, 

XIN( 

17)- 

0.4696 

,RIN( 

17)- 

7.1632, 

XIN( 

18)- 

0.5409 

,RIN( 

18)- 

7.1307, 

XIN( 

19)- 

0.6205 

,RIN( 

19)- 

7.0983, 

XIN( 

20)- 

0.7098 

,RIN( 

20)- 

7.0659, 


( Note - Complete 


listing of data not shown 


) 


XOUT( 82)- 35.0000 ,R0UT( 82)- 9.0000, 

X0UT( 83) -40. 0000 ,R0UT( 83)* 9.0000 $END 

$LIST3 $END 
$LIST4 $END 


Figure 30. Data deck for NGRIDA Sample Case No. 2. 
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listing 

OF 

INPUT DATA 


8LTST1 

IMAX 

3 

67. 


JMAX 

s 

13. 


KMAX 

s 

13. 


XLEFT 

s 

-30.00000 

. 

RAHCmT 

3 

50.00000 

. 

delta 

3 

O.bOOOOOU 


toutfr 

3 

38, 


tinner 

3 

30, 


alphai 

S 

i.iioooo 

9 

ALPHAO 

S 

1.125000 

9 

KOUTER 

3 

0, 


ALPHPI 

3 

1.000000 

9 

ALPHPO 

3 

1.000000 

9 

scale 

3 

U. 1190476 

9 

itape 

3 

10, 


KPLCJT 

3 

1. 


KDEV 

3 

1. 


KDUF.P 

3 

1. 


lAA 

3 

12. 


TDD 

3 

57. 


KPRIMT 

3 

1, 

24 

KTTMF 

3 

3. 


MDFL 

3 

5. 


DSFACT 

3 

1.000000 

9 

kramp 

3 

1. 


tramp 

3 

57 



SEND 

SLTST2 

NU = 8^, 

Xonj = O.OOOOOOOE+OO , 1 . 6 OOOOOOE-O 3 , 6 , 3000000E-03 , 1 . 4200000E-02 , 


2.53OOOOOE-02, 

3.9600000F-02, 

5.7100002F-02, 

7.7799998E-02, 

0.1017000 

$ 

0.1290000 , 

0.1596000 , 

0.1935000 

0.230800U , 

0.2716000 

$ 

0.3159000 , 

0.3638000 , 

0.4153000 

0.4705000 , 

0.5295000 

$ 


Figure 31. Output for NGRIDA Sample Case No. 2. 





CIRTD POINT COURUlfj ATKS pOk bP'LKUTtD MtRirin;jAl. PI.ANhS 


1 = 

1 

Ja 

1 

Ks 

1 

As 

O.H33i3 JF+01 

Ys 

0. bOS2 jPt; + 01 

Zs O.OOOOOOE-fOO 

Is 

2 

ja 

1 

Ks 

1 

As 

0.7rt5ildF+01 

Ys 

0 . 59b23flF + 01 

Zs n.noOoOOE-fOO 

Is 

3 

J* 

1 

Ks 

1 

Xs 

0.7 J«h29F+0l 

Ys 

0.b9b?3flF+u1 

Zs 0 , noOOOOE+OO 

Is 

4 

Ja 

1 

Ks 

1 

As 

0.69545b'^ + Ol 

Ys 

C . b9b?3Rfc>01 

Z= O.OoOOOOE^OO 

Is 


Ja 

1 

Ks 

1 

As 

O.6b670oF+01 

Ys 

0. b9b23PF + 01 

Zs O.OOOOOOE+OO 

Is 

6 

ja 

1 

Ks 

1 

As 

O.622129P' + 01 

Ys 

0, 59b?3Pt>01 

Zs O,O00000E+00 

Is 

7 

Ja 

1 

Ks 

1 

As 

0,5yl ORlF+01 

Ys 

L-, 59b23PE + 01 

Z= O.OOOOOOE+OO 

Is 

p 

ja 

1 

Ks 

1 

Xs 

0.5b2781E+01 

Ys 

0. b9b?38E+01 

Zs O.OoOOOOE+00 

Is 

9 

Ja 

1 

Ks 

1 

Xs 

0.5jftl22»=' + 01 

Ys 

C'.b952iRt>0l 

Zs O.OOOOOOE^OO 

Is 

10 

ja 

1 

Ks 

1 

Xa 

0, S0O41dr+01 

Ys 

C.b9b23RE-^0l 

Z= O.OOOOOOE^-OO 

Is 

11 

Ja 

1 

Ks 

1 

Xs 

0.4b0232t'+01 

Ys 

0.59b238fc,^01 

Zs O.OOOOOOE^OO 

Is 

12 

Ja 

1 

Ks 

1 

Xs 

0.445877E+01 

Ys 

0.59523PE+01 

Zs O.OOOoOOE+00 

Is 

13 

Ja 

1 

Ks 

1 

Xs 

0.405558F+01 

Ys 

C . b9b238E+0i 

Zs O.noOOOOE+OO 

Is 

14 

Ja 

1 

Ks 

1 

Xs 

0. 3b1 b7bF+0i 

Ys 

1 . £)9b?3Rt + 01 

Zs O.OOOOOOE+OO 

Is 

15 

Ja 

1 

Ks 

1 

Xs 

0. 3l6y9oF+01 

Ys 

0.b9b23Pt>01 

Z= O.ooOOOOE-fOO 

Is 

16 

J = 

1 

Ks 

1 

As 

0.272192F+01 

Ys 

C . 59b?38E+0l 

Zs O.OoOOOOE+OO 

Is 

17 

Js 

1 

Ks 

1 

As 

0.227077F+O1 

Ys 

C .bOb?38t + 0t 

Z= O.OOOOOOE+00 

Is 

IP 

Ja 

1 

Ks 

1 

As 

0. 17P964F+01 

Ys 

0. S9S23PE'^01 

Zs O.OOOOOOE-fOO 

Is 

19 

J = 

1 

Ks 

1 

Xs 

0. 127777F + 01 

Ys 

U,59b23RE+0t 

Z= O.OOOOOOE-fOO 

Is 

20 

J = 

1 

Ks 

1 

As 

0,7ei369F+C0 

Ys 

C. b87Q8PE*0l 

Zs O.OOOOOOE-fOO 

Is 

21 

J = 

1 

Ks 

1 

Xs 

0. 3d9639F+0U 

Ys 

( ,t>64275£ + 0l 

Z= O.ooOOOOE^OO 

Is 

22 

J = 

1 

Ks 

1 

As 

0.7b?49iF-01 

Ya 

0. bl b1 1 lE+OI 

Zs O.OOOOOOE+00 

Is 

23 

Ja 

1 

Ks 

1 

As 

0, OOOOOOK+00 

Ys 

0.45bOH3E+Ot 

Zs O.ooOOOOE-fOO 

Is 

24 

Ja 

1 

Ks 

1 

As 

O.OOOOOOF+00 

Ys 

C. 39b637F+01 

Zs O.OOOOOOE+OO 

I B 

25 

Ja 

1 

Ks 

1 

As 

O.noOOOOE+00 

Ys 

C.34b323E+01 

Zs O.noOOOOE+00 

Is 

26 

JS 

1 

Ks 

1 

As 

O.OOOOOOF+OO 

Ys 

C.301Q252+01 

Zs O.OOOOOOE+OO 

Is 

27 

J = 

1 

Ks 

1 

As 

O.OOOOOoFfOO 

Ys 

0.26lPd lE+01 

Z= o.noOOOOE+00 

Is 

2P 

Ja 

1 

Ks 

1 

As 

0, OUOOOOF^OO 

Ys 

C.27bO48E+01 

Zs O.ooOOOOE+00 

Is 

29 

Ja 

1 

Ks 

1 

As 

O.OOOOOOF-fOO 

Ys 

0, 1905b0E+01 

Z= O.OOOOOOE+OO 

T s 

30 

J = 

1 

Ks 

1 

As 

0,OOOOOOF+O0 

Ys 

0,15 /687£4.01 

Zs 0, OOOOOOE+OO 


Figure 31. Continued. 


omaiKkL PAGE m 

OF POOR QUALITY 



OF PWn QUALITY 


$LIST1 ALFH40-1.125, 
$LIST2 NI-64, NO-79, 

ALPHAI- 

1.11 

$END 

XIN( 

D- 

0.0000 

,RIH( 

D- 

7.6822, 

XIN( 

2)- 

0.0016 

,RIN( 

2)- 

7.6498, 

XIN( 

3)- 

0.0065 

,RIN( 

3)- 

7.6173, 

XIN( 

4)- 

0.0147 

,RIN( 

4)- 

7.5849, 

XIN( 

5)- 

0.0261 

,RIN( 

5)- 

7.5525, 

XIN( 

6)- 

0.0410 

,RIN( 

6)- 

7.5200, 

XIN( 

7)- 

0.0593 

,RIN( 

7)- 

7.4876, 

XIN( 

8)- 

0.0811 

,RIN( 

8)- 

7.4551, 

XIN( 

9)- 

0.1066 

,RIN( 

9)- 

7.4227, 

XIN( 

10)- 

0.1359 

,RIN( 

10)- 

7.3903, 

XIN( 

11)- 

0.1693 

,RIN( 

11)- 

7.3578, 

XIN( 

12)- 

0.2068 

,RIN( 

12)- 

7.3254, 

XIN( 

13)- 

0.2488 

,RIN( 

13)- 

7.2929, 

XIN( 

14)- 

0.2957 

,RIN( 

14)- 

7.2605, 

XIN( 

15)- 

0.3477 

,RIN( 

15)- 

7.2281, 

XIN( 

16)- 

0.4055 

,RIN( 

16)- 

7.1956, 

XIN( 

17)- 

0.4696 

,RIN( 

17)- 

7,1632, 

XIN( 

18)- 

0.5409 

,RIN( 

18)- 

7.1307, 

XIN( 

19)- 

0.6205 

,RIN( 

19)- 

7.0983, 

XIN( 

20)- 

0.7098 

,RIN( 

20)- 

7.0659, 


( Note - Complete listing of data not shown ) 


XOUT( 78)« 35.0000 ,R00T( 78)- 9.0000, 

X0UT( 79)- 40.0000 ,R0UT( 79)- 9.0000 $EKD 

$LIST3 $END 
$L1ST4 $END 


Figure 32. Data deck for NGRIDA Sample Case No. 3. 
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THRFfci-DIf.KwSION AL NACFLliE GRID GENERATION PROGRAM 


V_T\ 


LlSTTl^G 

OF 

INPUT DATA 

SLTST’l 



IMAX 

= 

67, 

JMAX 

s 

13, 

KMAX 

s 

13, 

xlfft 

s 

-3C.0O0O0 , 

RADOPT 

s 

50.00000 , 

DELIA 

s 

O.50C0OOO , 

lOUTFH 

s 

38. 

iinner 

s 

30. 

ALPHAI 

s 

I.IIOOOO , 

ALPHAO 

s 

1.125000 , 

KOUTFR 

s 

0. 

ALPhBI 

s 

l.OOOOOO , 

ALPHRU 

s 

1. 000000 , 

scale 

— 

0.1190476 , 

itapf 

s 

1^. 

KPLOT 

s 

1. 

kdev 

s 

1, 

KDPMP 

s 

t. 

lAA 

n 

12. 

IDD 

s 

57. 

KPRlMl 

3 

1. 24*0, 

KTTMF 

s 

3. 

NOEL 

s 

5, 

dsfact 

= 

1,000000 , 

KRAMP 

s 

1 . 

TRAMP 

s 

57 

$EMD 



8LIST2 



NO 

s 

79. 

XOl)T 

3 

o,oooooooF+no, 2 , 

3.1500001F 

-02, 4.9300000F-0 

0.1608000 

, 0.1988000 

0,3937000 

, 0,4533000 


OOOOOOlE-03, 7.BR99996F-0i, 1 . 7700000E-02 , 

2, 7, 1099997E-02, 9 . 6900v)0 1 E-02 , 0.1268000 

, 0,2411000 , 0,2877000 , 0.3385000 

, 0. 6175000 , 0,5863000 , 0.6599000 


* 

9 

9 


Figure 33. 


Output for NGRIDA Sample Case No. 3. 


OBIQSMAL PAGE si 
OF POOR QUALITY 



grid PniNX COURDlfJATEij FOR SFLECTtn MtRIDinNALi Pf.ANtS 


la 

1 

Ja 

1 

Ks 

1 

Xs 

0.833333F-f0l 

Ys 

0.595738fc:+01 

Zs 

o.oooooo£+oo 

la 

2 

Ja 

1 

Ka 

1 

Xa 

0,785363t>01 

Ys 

0.b9b?38fc:+01 

Za 

O.OOOOOOE+00 

la 

1 

ja 

1 

Ks 

1 

Xa 

0,738719F+ni 

Ys 

0.59b238t+01 

Za 

O.OOOOOOE+00 

la 

4 

J = 

1 

Ks 

1 

Xs. 

0.695585E-f01 

Ys 

a,59b?38E*01 

Za 

O.OOOOOOE+00 

la 

5 

Ja 

1 

Ks 

1 

Xs 

0.656873E+01 

Ys 

0.b9b238E+01 

Zs 

O.OOOOOOE+00 

la 

6 

JS 

1 

Ks 

1 

Xs 

0.622328E+01 

Ys 

0.b9b?38E+01 

Zs 

O.OOOOOOE+00 

T = 

7 

Ja 

1 

Ks 

1 

Xs 

0.591309E+01 

Ys 

0.59b?38E+01 

Za 

O.OOOOOOE+00 

la 

fl 

Ja 

1 

Ka 

1 

Xa 

0.563034F+01 

Ys 

0,b95238E+01 

Zs 

O.OOOOOOE+00 

Ta 

Q 

Ja 

1 

Ks 

1 

Xa 

0.53639oF+ni 

Ys 

0.b95738E+01 

Zs 

O.OOOOOOE+00 

la 

10 

Ja 

1 

Ks 

1 

Xs 

0.509707F+01 

Ys 

0,59b?38E+01 

Zs 

O.OOOOOOE+00 

la 

11 

Ja 

1 

Ks 

1 

Xs 

0.48052VF+01 

Ys 

0,b95?3RE+01 

Za 

O.OOOOOOE+00 

la 

12 

Ja 

1 

Ks 

1 

Xa 

0.446174E+01 

Ys 

0.59b73BE+01 

Z = 

O.OOOOOOE+00 

la 

n 

Ja 

1 

Ks 

1 

Xs 

0.405849F+01 

Ys 

0.b9b2J8E+01 

Za 

O.noOOOOE+00 

la 

14 

J* 

1 

Ka 

1 

Xa 

0,361958=' + 01 

Ys 

0.59573BE+01 

Za 

O.OOOOOOE+00 

la 

1«5 

Ja 

1 

Ks 

1 

Xa 

0.317272F+01 

Ys 

0.b95?38E+0t 

Zs 

O.OOOOOOE+00 

la 

If 

Ja 

1 

Ks 

1 

Xs 

0.27246bF+Ol 

Ys 

0.59b?3BE+01 

Za 

O.OOOOOOE+OO 

la 

17 

J* 

1 

Ks 

1 

Xs 

0,227350F+01 

Ys 

0,59b238E+01 

Zs 

O.OOOOOOE+00 

la 

18 

Ja 

1 

Ks 

1 

Xs 

0.17924bF+0l 

Ys 

0,b9b238E+01 

Za 

O.OOOOOOE+00 

la 

19 

J* 

1 

Ka 

1 

Xs 

0. 128071F + 01 

Ys 

0,b9b738E+01 

Za 

O.OOOOOOE+OO 

la 

20 

J = 

1 

Ks 

1 

Xs 

0.783893F+00 

Ys 

0.58bnBnE+01 

Zs 

O.OOOOOOE+00 

la 

21 

J = 

1 

Ka 

1 

Xs 

0.39186UF+no 

Ys 

U.b64477E+ut 

Za 

O.OOOOOOE+00 

Ta 

27 

JS 

1 

Ks 

t 

Xa 

0.77407ttF-01 

Ys 

O.bl 8418E+01 

Zs 

0 • OuOOOOE+00 

Is 

23 

J = 

1 

Ks 

1 

Xs 

0,OOOOOOF+O0 

Ya 

0,45b397E+Ul 

Za 

O.OOOOOOE+00 

Ta 

24 

JS 

1 

Ks 

1 

Xa 

O.OOOuOOF+OO 

Ys 

0.396P9SE+01 

Zs 

O.OOOOOOE+00 

la 

2*5 

JS 

1 

Ks 

1 

Xs 

O.OOOOOOF-fOO 

Ys 

0,3465<*1E + 01 

Za 

O.OOOOOOE+00 

Ta 

26 

JS 

1 

Ka 

1 

Xs 

O.OOOOOOF+00 

Ys 

0.302in£ + 01 

Zs 

O.OOOOOOE+OO 

la 

27 

JS 

1 

Ks 

1 

Xs 

o,oooooor+no 

Ya 

0.262O47E+0J 

Zs 

O.OOOOOOE+OO 

I a 

28 

Ja 

1 

Ks 

1 

Xs 

O.OOOOOOF+OO 

Ys 

0.225197E+01 

Z = 

O.OOOOOOE+00 

la 

29 

Js 

1 

Ks 

1 

Xs 

0, OOOOOOF+OO 

Ys 

0,190685E+01 

Zs 

O.OOOOOOE+OO 

la 

30 

JS 

1 

Ks 

J, 

Xs 

O.OOOOOOF+OO 

Ya 

0.157812E+01 

Zs 

0. noooooE +00 


Figure 33. Continued. 


OF POOR Q'UALiTV: 



PASS ^3- 

Of: POOR QUALITY 


$LIST1 ALFHAO-1.125, ALPHAI-l.U, JMAX-25, KMAX-16 $EHD 
$LIST2 RI-64, HO-83, 


XIN( 

D- 

0.0000 

,RIH( 

D- 

7.6822, 

XIN( 

2)- 

0.0016 

,RIH( 

2)- 

7.6498, 

XIN( 

3)- 

0.0065 

,RIN( 

3)- 

7.6173, 

XIN( 

4)- 

0.0147 

,RIN( 

4)- 

7.5849, 

XIN( 

5)- 

0.0261 

,RIN( 

5)- 

7.5525, 

XIN( 

6)- 

0.0410 

,RIN( 

6)- 

7.5200, 

XIN( 

7)- 

0.0593 

,RIN( 

7)- 

7.4876, 

XIN( 

8)- 

0.0811 

,RIN( 

8)- 

7.4551, 

XIN( 

9)- 

0.1066 

,RIN( 

9)- 

7.4227, 

XIN( 

10)- 

0.1359 

,RIN( 

10)- 

7.3903, 

XIN( 

11)- 

0.1693 

,RIN( 

11)- 

7.3578, 

XIN( 

12)- 

0.2068 

,RIN( 

12)- 

7.3254, 

XIN( 

13)- 

0.2488 

,RIN( 

13)- 

7.2929, 

XIN( 

14)- 

0.2957 

,RIN( 

14)- 

7.2605, 

XIN( 

15)- 

0.3477 

,RIN( 

15)- 

7.2281, 

XIN( 

16)» 

0.4055 

,RIN( 

16)- 

7.1956, 

XIN( 

17)- 

0.4696 

,RIN( 

17)- 

7.1632, 

XIN( 

18)- 

0.5409 

,RIN( 

18)- 

7.1307, 

XIN( 

19)- 

0.6205 

,RIN( 

19)- 

7.0983, 

XIN( 

20)- 

0.7098 

,RIN( 

20)- 

7.0659, 


( Note - Complete listing of data not shown ) 


XOUT( 82)- 35.0000 ,R0UT( 82)- 9.0000, 

X0UT( 83) -40. 0000 ,R0UT( 83)- 9.0000 $END 

$LIST3 $EHD 
$LIST4 $ERD 


Figure 34. Data deck for NGRIDA Sample Case No. 4. 
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OF POOR QUALITY, 

ALFHAO, ALFHAI, JMAZ» and KMAX. ALFHAO-I.125 and ALFHAI-I.IIO are again 
specified as before* In this case, however, the aesh resolution In the 
circumferential direction will be enhanced by entering JMAX«25 Instead of 
the default value JMAZ**13* Likewise the radial direction mesh resolution 
will be enhanced by specifying KMAX«16 Instead of the default value 
KMAX«13* The nacelle/ Inlet geometry data are again entered as before In 
KAMELIST LIST2. All parameters In NAMELISTS LIST3 and LIST4 are retained 
at their assigned default values* 

Selected portions of the computer output for Sample Case No. 4 are 
presented In Figure 35. 

3. NACELLE FLOW ANALYSIS PROGRAM SAMPLE CASES 

Four sample cases are presented in this section to illustrate 
application of the NACELLE flow solution algorithm to determining the flow 
field about realistic axisymmetric nacelle/ inlet configurations. The 
computational grid for each sample case is generated using the NGRIDA grid 
generation code. 

3.1 SAMPLE CASE NO. 1 

This sample case is concerned with determining the flow field solution 
for the Lockheed-Georgla GELACl axisymmetric nacelle/ inlet configuration at 
angle of attack. This sample case represents the case being considered 
when all input parameters in both the NACELLE flow solution algorithm and 
in the NGRIDA grid generation algorithm are retained at their default 
values. Consequently, the computational mesh for this sample case is 
generated by executing Sample Case No. 1 for the NGRIDA program. 

The data deck for Sample Case No. 1 is illustrated in Figure 36. All 
input parameters in NAMELIST LISTl are retained at their default values. 
Consequently, this case considers the flow computation for an angle of 
attack PITCH of 1.0 degree, a free-stream Mach number XMFS of 0.8, a 
compressor face Mach number XMCF of 0.3516, and a specific heat ratio GAMMA 
of 1.4. The computation will use second~order metrics and an internal flow 
Mach number initialization based on area ratio, as noted by the default 
values of K0RDER"0 and KINIT*!* The Internal flow Mach number multiplier 
XMULM is left at its default value of 0.7. The solution algorithm will 
stop when either the residual ratio reaches the value of 0.001 or when the 
iteration count reaches 200. The grid point coordinates will not be 
printed since KFRC0R*O, and the surface solution will be printed when 

the convergence criterion has been satisfied or when the iteration 
count reaches ITMAX, since KPRINT*1000. The printed surface coordinates 
will be up-scaled by the factor 8.4, since SCALE-8.4. This recaptures the 
original dimensions since scaling was performed in generating the grid. 

The Mach number and pressure coefficient distributions will be plotted for 
the J-1, J-4, and J-7 meridional planes, since KPLOT, KMACH, and KCP retain 
their default values. Since ITDUMP-0, the converged potential function 
field will not be loaded onto a disk file. 

All input parameters in NAMELIST LIST2 retain their default values. 
Consequently, the program assumes that a 67 x 13 x 13 (IMAX x JMAX x KMAX) 
computational mesh will be used. Since ITREAD-1 , the transition point 
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THRr.E-pIMKNSXDlM AL iiACFL^-t GPIH (iENERATIPN PROGRAM 


-n3 

NO 


listing 

OF 

INPUT DATA 

SLTSTl 

IMAX 

s 

67. 

JMAX 

B 

25, 

KMAX 

' s 

16, 

XLEFT 

B 

-30.00000 , 

RADOUT 

a 

50.00000 . 

delta 

a 

0,5000000 , 

rOUTFR 

s 

38, 

IIMNER 

a 

30, 

ALPHA! 

s 

1,110000 , 

ALPHAO 

a 

1.125000 , 

kouter 

a 

0, 

ALPHPI 

a 

l.OCOOOO , 

ALPHRO 

a 

1,000000 . 

scale 

a 

0.1190476 , 

itape 

a 

10, 

KPI.OT 

a 

1. 

KDEV 

a 

1 , 

koump 

s 

t. 

taa 

s 

17, 

IDO 

a 

57, 

KPRINT 

a 

1, 24*0, 

KTTMF 

a 

3, 

NDEL 

a 

5, 

nsFACT 

r 

1,000000 , 

kramp 

a 

1. 

IRAMP 

a 

57 

send 

SLTST2 

NO 

3 

83, 

XOUT 

3 

O.OOOOOOOE + OO , 1 

2,5100000E 

-02, 3.9600000E- 

0,1290000 

, 0,1596000 

0.3159000 

, 0.3638000 


lOOOE-0 3 , 6, 300000 OF. -03, 1 , 4 200Q00F-02 , 

5.7100002E-02, 7 , 7790998E-0;: , 0*10X7000 

0,1935000 , 0,2300000 , 0.2716000 

0,4153000 , 0.47U500O , 0.5295000 



Figure 35. Output for NGRIDA Sample Case No. 4. 


OF POOR Q 



00 

o 


Is 

3ft 

J = 

1 

Ka 

16 

Xs 0, 35761 dE+Ol 

Is 

37 

via 

1 

Ks 

Ift 

Xs 0.35725OF+01 

la 

3fl 

vJa 

1 

Ka 

Ift 

Xs 0.357143F: + Oi 

Is 

39 

Ja 

i 

Ks 

16 

Xa 0,357279F^01 

Is 

40 

Ja 
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Figure 35. Continued. 
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$LIST1 $END 
$LIST2 $END 
$LIST3 $END 
$LIST4 $END 


Figure 36. Data deck for NACELLE Sample Case No. 1. 
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indices ITRANl, ITRAN2, ITRAN3, and IBEG are read from the binary file 
produced by the grid generation algorithm. 

All Input parameters In NAMELIST LISTS retain their default values. 
Consequently, the sequence Is based on equation (39), ^nd 

o<^j-6.0. The relaxation factor OMEGA will be 1.0, and, since NjviS»NKVIS-0, 
upwlndlng will only be performed In the wraparound coordinate direction. 
CFACT and BXIE are left at their default values of 1.2 and 0.1, 
respectively. The smoothing function parameters KSMTH, KTIME, and NDEL are 
retained at their default values of 4, 3, and 3, respectively. The 
elements of the under-*relaxatlon factor array CORFAT are maintained at 
their monotonically increasing default values, starting at 0.2 and ending 
at I.O. 


No debug output is to be printed, hence all input parameters in 
NAMELIST LIST4 maintain their default values. 

Selected portions of the computer output for Sample Case No. 1 are 
presented in Figure 37. Results for the GSLACl nacelle are presented in 
Section II. 

3.2 SAMPLE CASE NO. 2 

This sample case is concerned with determining the flow field solution 
for the Lockheed-Californla CALAC4 axisymmetric nacelle/ inlet configuration 
at incidence. The computational mesh for this case is generated by 
executing Sample Case No. 2 for the NGRIDA program. 

The data deck for this sample case is illustrated in Figure 38. All 
input parameters in NAMELISTS LIST2 and LIST4 retain their default values. 
In this case the angle of attack is 2.0 degrees, hence PITCH»2.0 is 
specified in NAMELIST LISTl. The remaining parameters in NAMELIST LISTl 
are left at their default values. All parameters in NAMELIST LIST3 retain 
their default values except for CFACT, which has now been slightly 
increased to 1.3. 

Selected portions of the computer output for Sample Case No. 2 are 
presented in Figure 39. Results for the CALAC4 nacelle are presented in 
Section II. 

3.3 SAMPLE CASE NO. 3 

This sample case is concerned with the flow field computation for the 
Lockheed-California CALAC3 axisymmetric nacelle geometry. This solution is 
for zero incidence, hence the flow field will be axisymmetric. The grid 
for this case is determined by executing Sample Case No. 3 for the NGRIDA 
program. 

The data deck for this sample case is illustrated in Figure 40. All 
input parameters retain their default values except for PITCH in NAMELIST 
LISTl. PITCH is now equated to 0.0. 

Selected portions of the computer output for Sample Case No. 3 are 
presented in Figure 41. 
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Figure 37. 

Output for NACELLE Sample Case No. 1. 
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0.7130 

MACHs 

0,851 1 

cps -0.1130 

-n ^ 

Is 

1^ 

JS 

X 

Xs 

3ft .0746 

nxs 

d.874V> 

Ya 

9.00CO 

z= 

0, onoo 

RHJs 

0.70b 1 

MACHs 

0 . ftb03 

rps -0.1334 

-e 0 

1 = 

Ift 

JS 

X 

Xs 

J7.Q243 

nxs 

7.9245 

Ya 

a ,9095 

Zs 

0,0000 

RhU- 

0,6999 

y.ACHs 

0.8756 

CPs -0.1674 

0 W 

Is 

17 

JS 

X 

Xs 

37 , ObOd 

nxs 

7.0600 

Ya 

a . 9-500 

z= 

O.OOuD 

RH0 = 

0.6877 

•SACHS 
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$LIST1 PITCH=2.0 $END 
$LIST2 $END 
$LIST3 CFACT=1.3 $END 
$LIST4 $END 


Figure 38. Data deck for NACELLE Sample Case No. 2. 
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LISTING OF INPUT DATA 
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Figure 39. Output for NACELLE Sample Case No, 2. 
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SURFACF SOLUTION FDR ITERATinN NO.® 166 
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Figure 40. Data deck for NACELLE Sample Case No. 3. 
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Figure 41. Output for NACELLE Sample Case No. 3. 
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3.4 SAMPLE CASE NO. 4 
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This sample case is concerned with the flow field computation for the 
CALAC4 nacelle using a (67 X 25 X 16)-point grid. The computational mesh 
for this case is generated by executing Sample Case No. 4 for the NGRXDA 
program. 

The data deck for this sample case is illustrated in Figure 42. All 
input parameters in NAMELIST LISTl retain their default values except for 
XMFS, CRATIO, PITCH, ITMAX, KMACH, and KCP. For this case, the free“stream 
Mach number XMFS is 0.6, the capture ratio CRATIO is 0.505, and the angle 
of attack PITCH is 1.083 degrees. To ensure convergence, ITMAX has been 
increased to 300. The surface solution is to be plotted for the meridians 
at 0.0 degrees, 90.0 degrees, and 180.0 degrees; hence, 

KMACH(1)=1, 5*0, 1,5*0, 1,12*0 and KCP(1 )=1 , 5*0 , 1 , 5*0 , 1 , 12*0 are entered in 
NAliELIST LISTl. 

Since a (67 X 25 X 16)-point grid will be employed, JMAX=25 and KMAX=16 
are entered in NAMELIST LIST2. IMAX is retained at its default value of 
67. 


All input parameters in NAMELIST LIST3 retain their default values 
except for ALPHAL, ALPHAH, and CORF AT. To optimize convergence speed, 
ALPHAL0.4 and ALPHAH=4.0 are entered. These values were obtained by 
numerical experiment. Since 16 radial stations are to be used, the 
elements of the CORFAT array need to be altered from their default values. 
This array is now specified as C0RFAT(K)=. 2, .25, .3, .35, .4, .45, .5, .6, 
.7, .8, .9, and 5*1.0 for K=1 to K=16, respectively. Note that 
C0RFAT(1)=.2 and CORFAT (KMAX)=*1 .0 as in the default assignment. 

All parameters in NAMELIST LIST4 retain their default assignments since 
no debug output is to be printed. 

Selected portions for the computer output for Sample Case No. 4 are 
presented in Figure 43. Results for the CALAC4 nacelle are presented in 
Section II. 

4. SUGGESTIONS FOR USAGE 
4.1 PARAMETER STATEMENT USAGE 

The NGRIDA grid generation program and the NACELLE flow analysis 
program have been executed on the VAX-1 1/7 80 and CRAY-1 computers. Both 
the NGRIDA and NACELLE programs employ PARAMETER statements which fix the 
respective array dimensions at the time of compilation. This provides an 
effective means of changing program core storage requirements. The VAX 
version of each program employs an INCLUDE statement which inserts the 
respective PARAMETER statement at the appropriate locations within the 
code. The NGRIDA and NACELLE programs employ the ’INCLUDE 'NGRID.DIM' ' and 
’INCLUDE ’NACELLE . DIM ’ ’ statements, respectively. If the codes are 
executed on a VAX, then the NGRID.DIM and NACELLE. DIM files must reside in 
the user’s directory. For execution on the CRAY-1 or CDC computers, the 
PARAMETER statement is edited into the code by use of a COMDECK. For use 
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$LIST1 XMFS-0.6, CRATIO-0.505, PITCH=1.083, ITMAX*=300, 

KMACH(1)-1, 5*0, 1, 5*0, 1, 12*0, 

KCP(1)-1, 5*0, 1, 5*0, 1, 12*0 $END 
$LIST2 KMAX«16, JMAX-25 $END 
$LIST3 ALPHAL-0.4, ALPHAH-4.0, 

CORFAT(l)- .2, .25, .3, .35, ,4, .45, .5, ,6, .7, ,8, .9, 5*1.0 $END 
$LIST4 $END 


Figure 42. Data deck for NACELLE Sample Case No. 4. 
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OF INPUT DATA 
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IbFG 

s 

0 


SEND 

SLIST3 

KALPHA 


t . 


ALPHAL 

S 

O.400O00O 

0 

AlPhAH 

- 
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Figure 43. Output for NACELLE Sample Case No. 4 
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on other conputer systeas, the PARAMETER etatenents will have to be 
inserted in the respective programs at the locations noted by the INCLUDE 
statenents* 

4.2 NOLIDA PROQLAM USAGE 

The NGRIDA grid generation program often can be executed by retaining 
most of the input parameters at their respective default values. 

Maintaining the outer computational boundary at the same distance from the 
body, relative to the compressor face radius (RCF), as was done In the 
supplied sample cases should be sufficient for most new applications of the 
program. This can be accomplished by maintaining the ratios RADOUT/RCF 
XLEFT/RCF, and DELTA/RCF the same as in the sample cases provided herein. 

For transonic flow computations, supercritical flow occurs generally in 
the external flow region that is slightly downstream of the inlet lip. The 
mesh resolution can be enhanced in this region by increasing the mesh 
stretching parameter ALPHAO beyond its default value of L.l. The C-mesh 
spacing is quite sensitive to ALPHAO and only a small Increase beyond its 
default value is usually required. Increasing the number of 
circumferential stations, denoted by JMAX, generally improves the ability 
of the flow analysis program to predict the Mach number distribution peaks 
especially on the external surface of the leeward meridian. 

In the sample cases provided herein, the grid coordinates are 
down-scaled by setting the SCALE input parameter equal to the reciprocal of 
the compressor face radius. This practice is recommended as it produces a 
computational mesh with a compressor radius of unity and thereby sizes the 
grid so that the default AF2 algorithm acceleration parameters, denoted by 
ALPHA! and ALPHAH in the NACELLE program input , are approximately optimum. 
The original coordinates can then be recovered in the NACELLE flow solution 
program by setting the SCALE input parameter in that code equal to the 
compressor face radius. 


T grape algorithm input parameters specified in NAMELISTS 

LISTS and LIST4 of the NGRIDA program are retained at their default values. 
If a nacelle with a very sharp inlet lip is to be analyzed, however, it may 
be necessary to reduce or equate to 0.0 the values of OMEGP and OMEGO in 
NAMELIST LIST4. This relaxes the enforcement of the grid spacing and 
orthogonality conditions at the nacelle surface and should permit a 
converged grid solution to be obtained. 

4.3 NACELLE PROGRAM USAGE 


The NACELLE flow analysis program often can be executed by retaining 
most of the input parameters at their respective default values. For cases 
with very strong shocks or high local supersonic Mach numbers, the 
artificial viscosity parameter CFACT and the time-like dissipation factor 
BXIE will probably have to be Increased over their default values of 1.2 
and 0.1, respectively. Increasing CFACT while holding BXIE at its default 
value for a typical case can decrease the number of Iterations required for 
convergence. However, increasing CFACT generally decreases shock wave 
resolution which is indicated by a slight smearing of the surface pressure 
distribution in the vicinity of the shock. Increasing BXIE generally 
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decreases convergence speed but ensures numerical stability by maintaining 
diagonal domi n a n ce In the ^-*factorlzed difference equations. 

The AF2 acceleration parameters ALFHAL and ALPHAH and the potential 
correction under--relaxatlon array CORFAT will require Input values 
different from the default values If different size grids are employed. 

The default values of ALFHAm).l75 and ALPHAH*^. 0 were obtained by 
numerical experiment to produce the optimum convergence speed on the 
(67 X 13 X 13)~polnt default computational mesh. If finer meshes are to be 
used, ALPHAL and ALPHAH will have to be altered. For a (67 x 25 x 
16)*‘polnt grid, ALFHAL*0.4 and ALPHAH«4.0 provided good convergence speed. 
Because of the non-linear nature of the analysis. It Is not possible to 
predetermine the optimum values of the ot parameters, hence they generally 
are determined by numerical experiment for a given grid size. Minor 
alterations in the surface geometry generally have little effect on 
convergence speed once the acceleration parameters are optimized for a 
given grid size. 

If the number of radial stations, denoted by KMAX, is changed then so 
must the CORFAT potential correction under-relaxation array. It is 
recommended that CORFAT(l) ~ .15, .2 and that CORFAT (KMAX )»1 . 0. The 
remaining elements of the CORFAT array should vary between these extremes 
and should monotonically increase with increasing K-lndex. Generally, 
fastest convergence is attained if the maximum residual occurs on the body 
surface and not near the centerline. If it occurs near the centerline, 
then slightly increasing the first few elements of the CORFAT array should 
move it to the body. 
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